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

Last change on this file since 37 was 37, checked in by raasch, 15 years ago

new version number 3.1c, prandtl layer switch used in production_e

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