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

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

Initial repository layout and content

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