Home


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

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

Last change on this file since 109 was 108, checked in by letzel, 17 years ago
  • Improved coupler: evaporation - salinity-flux coupling for humidity = .T.,

avoid MPI hangs when coupled runs terminate, add DOC/app/chapter_3.8;

  • Optional calculation of km and kh from initial TKE e_init;
  • Default initialization of km,kh = 0.00001 for ocean = .T.;
  • Allow data_output_pr= q, wq, w"q", w*q* for humidity = .T.;
  • Bugfix: Rayleigh damping for ocean fixed.
  • Property svn:keywords set to Id
File size: 40.5 KB
Line 
1 MODULE production_e_mod
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6! Bugfix: wrong sign removed from the buoyancy production term in the case
7! use_reference = .T.,
8! u_0 and v_0 are calculated for nxr+1, nyn+1 also (otherwise these values are
9! not available in case of non-cyclic boundary conditions)
10! Bugfix for ocean density flux at bottom
11!
12! Former revisions:
13! -----------------
14! $Id: production_e.f90 108 2007-08-24 15:10:38Z letzel $
15!
16! 97 2007-06-21 08:23:15Z raasch
17! Energy production by density flux (in ocean) added
18! use_pt_reference renamed use_reference
19!
20! 75 2007-03-22 09:54:05Z raasch
21! Wall functions now include diabatic conditions, call of routine wall_fluxes_e,
22! reference temperature pt_reference can be used in buoyancy term,
23! moisture renamed humidity
24!
25! 37 2007-03-01 08:33:54Z raasch
26! Calculation extended for gridpoint nzt, extended for given temperature /
27! humidity fluxes at the top, wall-part is now executed in case that a
28! Prandtl-layer is switched on (instead of surfaces fluxes switched on)
29!
30! RCS Log replace by Id keyword, revision history cleaned up
31!
32! Revision 1.21  2006/04/26 12:45:35  raasch
33! OpenMP parallelization of production_e_init
34!
35! Revision 1.1  1997/09/19 07:45:35  raasch
36! Initial revision
37!
38!
39! Description:
40! ------------
41! Production terms (shear + buoyancy) of the TKE
42! WARNING: the case with prandtl_layer = F and use_surface_fluxes = T is
43!          not considered well!
44!------------------------------------------------------------------------------!
45
46    USE wall_fluxes_mod
47
48    PRIVATE
49    PUBLIC production_e, production_e_init
50
51    LOGICAL, SAVE ::  first_call = .TRUE.
52
53    REAL, DIMENSION(:,:), ALLOCATABLE, SAVE ::  u_0, v_0
54
55    INTERFACE production_e
56       MODULE PROCEDURE production_e
57       MODULE PROCEDURE production_e_ij
58    END INTERFACE production_e
59   
60    INTERFACE production_e_init
61       MODULE PROCEDURE production_e_init
62    END INTERFACE production_e_init
63 
64 CONTAINS
65
66
67!------------------------------------------------------------------------------!
68! Call for all grid points
69!------------------------------------------------------------------------------!
70    SUBROUTINE production_e
71
72       USE arrays_3d
73       USE cloud_parameters
74       USE control_parameters
75       USE grid_variables
76       USE indices
77       USE statistics
78
79       IMPLICIT NONE
80
81       INTEGER ::  i, j, k
82
83       REAL    ::  def, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, &
84                   k1, k2, theta, temp
85
86!       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  usvs, vsus, wsus, wsvs
87       REAL, DIMENSION(nzb:nzt+1) ::   usvs, vsus, wsus, wsvs
88
89!
90!--    First calculate horizontal momentum flux u'v', w'v', v'u', w'u' at
91!--    vertical walls, if neccessary
92!--    So far, results are slightly different from the ij-Version.
93!--    Therefore, ij-Version is called further below within the ij-loops.
94!       IF ( topography /= 'flat' )  THEN
95!          CALL wall_fluxes_e( usvs, 1.0, 0.0, 0.0, 0.0, wall_e_y )
96!          CALL wall_fluxes_e( wsvs, 0.0, 0.0, 1.0, 0.0, wall_e_y )
97!          CALL wall_fluxes_e( vsus, 0.0, 1.0, 0.0, 0.0, wall_e_x )
98!          CALL wall_fluxes_e( wsus, 0.0, 0.0, 0.0, 1.0, wall_e_x )
99!       ENDIF
100
101!
102!--    Calculate TKE production by shear
103       DO  i = nxl, nxr
104
105          DO  j = nys, nyn
106             DO  k = nzb_diff_s_outer(j,i), nzt
107
108                dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
109                dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
110                                 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
111                dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
112                                 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
113
114                dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
115                                 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
116                dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
117                dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
118                                 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
119
120                dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
121                                 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
122                dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
123                                 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
124                dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
125
126                def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
127                      dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
128                      dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
129
130                IF ( def < 0.0 )  def = 0.0
131
132                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
133               
134             ENDDO
135          ENDDO
136
137          IF ( prandtl_layer )  THEN
138
139!
140!--          Position beneath wall
141!--          (2) - Will allways be executed.
142!--          'bottom and wall: use u_0,v_0 and wall functions'
143             DO  j = nys, nyn
144
145                IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) &
146                THEN
147
148                   k = nzb_diff_s_inner(j,i) - 1
149                   dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
150                   dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
151                                  u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
152                   dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
153                   dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
154                                  v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
155                   dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
156
157                   IF ( wall_e_y(j,i) /= 0.0 )  THEN
158                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
159                                          usvs, 1.0, 0.0, 0.0, 0.0 )
160                      dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
161!                      dudy = wall_e_y(j,i) * usvs(k,j,i) / km(k,j,i)
162                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
163                                          wsvs, 0.0, 0.0, 1.0, 0.0 )
164                      dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
165!                      dwdy = wall_e_y(j,i) * wsvs(k,j,i) / km(k,j,i)
166                   ELSE
167                      dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
168                                      u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
169                      dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
170                                      w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
171                   ENDIF
172
173                   IF ( wall_e_x(j,i) /= 0.0 )  THEN
174                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
175                                          vsus, 0.0, 1.0, 0.0, 0.0 )
176                      dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
177!                      dvdx = wall_e_x(j,i) * vsus(k,j,i) / km(k,j,i)
178                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
179                                          wsus, 0.0, 0.0, 0.0, 1.0 )
180                      dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
181!                      dwdx = wall_e_x(j,i) * wsus(k,j,i) / km(k,j,i)
182                   ELSE
183                      dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
184                                      v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
185                      dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
186                                      w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
187                   ENDIF
188
189                   def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
190                         dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
191                         dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
192
193                   IF ( def < 0.0 )  def = 0.0
194
195                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
196
197
198!
199!--                (3) - will be executed only, if there is at least one level
200!--                between (2) and (4), i.e. the topography must have a
201!--                minimum height of 2 dz. Wall fluxes for this case have
202!--                already been calculated for (2).
203!--                'wall only: use wall functions'
204
205                   DO  k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2
206
207                      dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
208                      dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
209                                     u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
210                      dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
211                      dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
212                                     v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
213                      dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
214
215                      IF ( wall_e_y(j,i) /= 0.0 )  THEN
216                         dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
217!                         dudy = wall_e_y(j,i) * usvs(k,j,i) / km(k,j,i)
218                         dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
219!                         dwdy = wall_e_y(j,i) * wsvs(k,j,i) / km(k,j,i)
220                      ELSE
221                         dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
222                                         u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
223                         dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
224                                         w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
225                      ENDIF
226
227                      IF ( wall_e_x(j,i) /= 0.0 )  THEN
228                         dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
229!                         dvdx = wall_e_x(j,i) * vsus(k,j,i) / km(k,j,i)
230                         dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
231!                         dwdx = wall_e_x(j,i) * wsus(k,j,i) / km(k,j,i)
232                      ELSE
233                         dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
234                                         v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
235                         dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
236                                         w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
237                      ENDIF
238
239                      def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
240                           dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 +  &
241                           dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
242
243                      IF ( def < 0.0 )  def = 0.0
244
245                      tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
246
247                   ENDDO
248
249                ENDIF
250
251             ENDDO
252
253!
254!--          (4) - will allways be executed.
255!--          'special case: free atmosphere' (as for case (0))
256             DO  j = nys, nyn
257
258                IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) &
259                THEN
260
261                   k = nzb_diff_s_outer(j,i)-1
262
263                   dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
264                   dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
265                                    u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
266                   dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
267                                    u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
268
269                   dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
270                                    v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
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
275                   dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
276                                    w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
277                   dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
278                                    w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
279                   dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
280
281                   def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
282                         dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
283                         dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
284
285                   IF ( def < 0.0 )  def = 0.0
286
287                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
288
289                ENDIF
290
291             ENDDO
292
293!
294!--          Position without adjacent wall
295!--          (1) - will allways be executed.
296!--          'bottom only: use u_0,v_0'
297             DO  j = nys, nyn
298
299                IF ( ( wall_e_x(j,i) == 0.0 ) .AND. ( wall_e_y(j,i) == 0.0 ) ) &
300                THEN
301
302                   k = nzb_diff_s_inner(j,i)-1
303
304                   dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
305                   dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
306                                    u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
307                   dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
308                                    u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
309
310                   dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
311                                    v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
312                   dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
313                   dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
314                                    v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
315
316                   dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
317                                    w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
318                   dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
319                                    w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
320                   dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
321
322                   def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
323                         dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
324                         dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
325
326                   IF ( def < 0.0 )  def = 0.0
327
328                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
329               
330                ENDIF
331
332             ENDDO
333
334          ELSEIF ( use_surface_fluxes )  THEN
335
336             DO  j = nys, nyn
337
338                k = nzb_diff_s_outer(j,i)-1
339
340                dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
341                dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
342                                 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
343                dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
344                                 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
345
346                dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
347                                 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
348                dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
349                dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
350                                 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
351
352                dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
353                                 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
354                dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
355                                 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
356                dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
357
358                def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
359                      dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
360                      dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
361
362                IF ( def < 0.0 )  def = 0.0
363
364                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
365
366             ENDDO
367
368          ENDIF
369
370!
371!--       Calculate TKE production by buoyancy
372          IF ( .NOT. humidity )  THEN
373
374             IF ( use_reference )  THEN
375
376                IF ( ocean )  THEN
377!
378!--                So far in the ocean no special treatment of density flux in
379!--                the bottom and top surface layer
380                   DO  j = nys, nyn
381                      DO  k = nzb_s_inner(j,i), nzt
382                         tend(k,j,i) = tend(k,j,i) +                    &
383                                       kh(k,j,i) * g / prho_reference * &
384                                       ( rho(k+1,j,i)-rho(k-1,j,i) ) * dd2zu(k)
385                      ENDDO
386                   ENDDO
387
388                ELSE
389
390                   DO  j = nys, nyn
391                      DO  k = nzb_diff_s_inner(j,i), nzt_diff
392                         tend(k,j,i) = tend(k,j,i) -                  &
393                                       kh(k,j,i) * g / pt_reference * &
394                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
395                      ENDDO
396
397                      IF ( use_surface_fluxes )  THEN
398                         k = nzb_diff_s_inner(j,i)-1
399                         tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
400                      ENDIF
401
402                      IF ( use_top_fluxes )  THEN
403                         k = nzt
404                         tend(k,j,i) = tend(k,j,i) + g / pt_reference * &
405                                                     tswst(j,i)
406                      ENDIF
407                   ENDDO
408
409                ENDIF
410
411             ELSE
412
413                IF ( ocean )  THEN
414!
415!--                So far in the ocean no special treatment of density flux in
416!--                the bottom and top surface layer
417                   DO  j = nys, nyn
418                      DO  k = nzb_s_inner(j,i), nzt
419                         tend(k,j,i) = tend(k,j,i) -                &
420                                       kh(k,j,i) * g / rho(k,j,i) * &
421                                       ( rho(k+1,j,i)-rho(k-1,j,i) ) * dd2zu(k)
422                      ENDDO
423                   ENDDO
424
425                ELSE
426
427                   DO  j = nys, nyn
428                      DO  k = nzb_diff_s_inner(j,i), nzt_diff
429                         tend(k,j,i) = tend(k,j,i) -               &
430                                       kh(k,j,i) * g / pt(k,j,i) * &
431                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
432                      ENDDO
433
434                      IF ( use_surface_fluxes )  THEN
435                         k = nzb_diff_s_inner(j,i)-1
436                         tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
437                      ENDIF
438
439                      IF ( use_top_fluxes )  THEN
440                         k = nzt
441                         tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
442                      ENDIF
443                   ENDDO
444
445                ENDIF
446
447             ENDIF
448
449          ELSE
450
451             DO  j = nys, nyn
452
453                DO  k = nzb_diff_s_inner(j,i), nzt_diff
454
455                   IF ( .NOT. cloud_physics )  THEN
456                      k1 = 1.0 + 0.61 * q(k,j,i)
457                      k2 = 0.61 * pt(k,j,i)
458                   ELSE
459                      IF ( ql(k,j,i) == 0.0 )  THEN
460                         k1 = 1.0 + 0.61 * q(k,j,i)
461                         k2 = 0.61 * pt(k,j,i)
462                      ELSE
463                         theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
464                         temp  = theta * t_d_pt(k)
465                         k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
466                                    ( q(k,j,i) - ql(k,j,i) ) *          &
467                              ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
468                              ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
469                              ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
470                         k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
471                      ENDIF
472                   ENDIF
473
474                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * &
475                                      ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
476                                        k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
477                                      ) * dd2zu(k)
478                ENDDO
479
480             ENDDO
481
482             IF ( use_surface_fluxes )  THEN
483
484                DO  j = nys, nyn
485
486                   k = nzb_diff_s_inner(j,i)-1
487
488                   IF ( .NOT. cloud_physics )  THEN
489                      k1 = 1.0 + 0.61 * q(k,j,i)
490                      k2 = 0.61 * pt(k,j,i)
491                   ELSE
492                      IF ( ql(k,j,i) == 0.0 )  THEN
493                         k1 = 1.0 + 0.61 * q(k,j,i)
494                         k2 = 0.61 * pt(k,j,i)
495                      ELSE
496                         theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
497                         temp  = theta * t_d_pt(k)
498                         k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
499                                    ( q(k,j,i) - ql(k,j,i) ) *          &
500                              ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
501                              ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
502                              ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
503                         k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
504                      ENDIF
505                   ENDIF
506
507                   tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
508                                         ( k1* shf(j,i) + k2 * qsws(j,i) )
509                ENDDO
510
511             ENDIF
512
513             IF ( use_top_fluxes )  THEN
514
515                DO  j = nys, nyn
516
517                   k = nzt
518
519                   IF ( .NOT. cloud_physics )  THEN
520                      k1 = 1.0 + 0.61 * q(k,j,i)
521                      k2 = 0.61 * pt(k,j,i)
522                   ELSE
523                      IF ( ql(k,j,i) == 0.0 )  THEN
524                         k1 = 1.0 + 0.61 * q(k,j,i)
525                         k2 = 0.61 * pt(k,j,i)
526                      ELSE
527                         theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
528                         temp  = theta * t_d_pt(k)
529                         k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
530                                    ( q(k,j,i) - ql(k,j,i) ) *          &
531                              ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
532                              ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
533                              ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
534                         k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
535                      ENDIF
536                   ENDIF
537
538                   tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
539                                         ( k1* tswst(j,i) + k2 * qswst(j,i) )
540                ENDDO
541
542             ENDIF
543
544          ENDIF
545
546       ENDDO
547
548    END SUBROUTINE production_e
549
550
551!------------------------------------------------------------------------------!
552! Call for grid point i,j
553!------------------------------------------------------------------------------!
554    SUBROUTINE production_e_ij( i, j )
555
556       USE arrays_3d
557       USE cloud_parameters
558       USE control_parameters
559       USE grid_variables
560       USE indices
561       USE statistics
562
563       IMPLICIT NONE
564
565       INTEGER ::  i, j, k
566
567       REAL    ::  def, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, &
568                   k1, k2, theta, temp
569
570       REAL, DIMENSION(nzb:nzt+1) ::  usvs, vsus, wsus, wsvs
571
572!
573!--    Calculate TKE production by shear
574       DO  k = nzb_diff_s_outer(j,i), nzt
575
576          dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
577          dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
578                           u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
579          dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
580                           u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
581
582          dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
583                           v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
584          dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
585          dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
586                           v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
587
588          dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
589                           w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
590          dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
591                           w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
592          dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
593
594          def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 )                       &
595                + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 &
596                + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
597
598          IF ( def < 0.0 )  def = 0.0
599
600          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
601               
602       ENDDO
603
604       IF ( prandtl_layer )  THEN
605
606          IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) )  THEN
607
608!
609!--          Position beneath wall
610!--          (2) - Will allways be executed.
611!--          'bottom and wall: use u_0,v_0 and wall functions'
612             k = nzb_diff_s_inner(j,i)-1
613
614             dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
615             dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
616                            u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
617             dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
618             dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
619                            v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
620             dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
621
622             IF ( wall_e_y(j,i) /= 0.0 )  THEN
623                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
624                                    usvs, 1.0, 0.0, 0.0, 0.0 )
625                dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
626                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
627                                    wsvs, 0.0, 0.0, 1.0, 0.0 )
628                dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
629             ELSE
630                dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
631                                u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
632                dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
633                                w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
634             ENDIF
635
636             IF ( wall_e_x(j,i) /= 0.0 )  THEN
637                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
638                                    vsus, 0.0, 1.0, 0.0, 0.0 )
639                dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
640                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
641                                    wsus, 0.0, 0.0, 0.0, 1.0 )
642                dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
643             ELSE
644                dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
645                                v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
646                dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
647                                w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
648             ENDIF
649
650             def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
651                   dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
652                   dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
653
654             IF ( def < 0.0 )  def = 0.0
655
656             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
657
658!
659!--          (3) - will be executed only, if there is at least one level
660!--          between (2) and (4), i.e. the topography must have a
661!--          minimum height of 2 dz. Wall fluxes for this case have
662!--          already been calculated for (2).
663!--          'wall only: use wall functions'
664             DO  k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2
665
666                dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
667                dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
668                               u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
669                dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
670                dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
671                               v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
672                dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
673
674                IF ( wall_e_y(j,i) /= 0.0 )  THEN
675                   dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
676                   dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
677                ELSE
678                   dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
679                                   u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
680                   dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
681                                   w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
682                ENDIF
683
684                IF ( wall_e_x(j,i) /= 0.0 )  THEN
685                   dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
686                   dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
687                ELSE
688                   dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
689                                   v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
690                   dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
691                                   w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
692                ENDIF
693
694                def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
695                      dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
696                      dvdz**2 + 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!
705!--          (4) - will allways be executed.
706!--          'special case: free atmosphere' (as for case (0))
707             k = nzb_diff_s_outer(j,i)-1
708
709             dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
710             dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
711                              u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
712             dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
713                              u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
714
715             dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
716                              v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
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(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
720
721             dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
722                              w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
723             dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
724                              w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
725             dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
726
727             def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
728                   dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
729                   dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
730
731             IF ( def < 0.0 )  def = 0.0
732
733             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
734
735          ELSE
736
737!
738!--          Position without adjacent wall
739!--          (1) - will allways be executed.
740!--          'bottom only: use u_0,v_0'
741             k = nzb_diff_s_inner(j,i)-1
742
743             dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
744             dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
745                              u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
746             dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
747                              u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
748
749             dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
750                              v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
751             dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
752             dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
753                              v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
754
755             dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
756                              w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
757             dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
758                              w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
759             dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
760
761             def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 )                       &
762                   + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 &
763                   + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
764
765             IF ( def < 0.0 )  def = 0.0
766
767             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
768
769          ENDIF
770
771       ELSEIF ( use_surface_fluxes )  THEN
772
773          k = nzb_diff_s_outer(j,i)-1
774
775          dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
776          dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
777                           u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
778          dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
779                           u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
780
781          dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
782                           v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
783          dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
784          dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
785                           v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
786
787          dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
788                           w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
789          dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
790                           w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
791          dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
792
793          def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
794                dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
795                dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
796
797          IF ( def < 0.0 )  def = 0.0
798
799          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
800
801       ENDIF
802
803!
804!--    Calculate TKE production by buoyancy
805       IF ( .NOT. humidity )  THEN
806
807          IF ( use_reference )  THEN
808
809             IF ( ocean )  THEN
810!
811!--             So far in the ocean no special treatment of density flux in the
812!--             bottom and top surface layer
813                k= nzb_s_inner(j,i)
814                tend(k,j,i) = tend(k,j,i) + kh(k,j,i) * g / prho_reference * &
815                     ( rho(k+1,j,i) - rho(k,j,i) ) * dzu(k+1)
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), 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.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |