source: palm/tags/release-3.1c/SOURCE/production_e.f90 @ 43

Last change on this file since 43 was 39, checked in by raasch, 17 years ago

comments prepared for 3.1c

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