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

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

preliminary version of modified boundary conditions at top

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