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

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

preliminary version, several changes to be explained later

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