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

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

RCS Log replace by Id keyword, revision history cleaned up

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