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

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

small bugs removed

  • Property svn:keywords set to Id
File size: 34.7 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 55 2007-03-08 04:12:41Z 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 beneath wall
114!--          (2) - Will allways be executed.
115!--          'bottom and wall: use u_0,v_0 and 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) - will be executed only, if there is at least one level
170!--                between (2) and (4), i.e. the topography must have a
171!--                minimum height of 2 dz. Wall fluxes for this case have
172!--                already been calculated for (2).
173!--                'wall only: use wall functions'
174
175                   DO  k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2
176
177                      dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
178                      dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
179                                     u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
180                      dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
181                      dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
182                                     v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
183                      dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
184
185                      IF ( wall_e_y(j,i) /= 0.0 )  THEN
186                         dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
187                         dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
188                      ELSE
189                         dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
190                                         u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
191                         dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
192                                         w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
193                      ENDIF
194
195                      IF ( wall_e_x(j,i) /= 0.0 )  THEN
196                         dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
197                         dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
198                      ELSE
199                         dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
200                                         v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
201                         dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
202                                         w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
203                      ENDIF
204
205                      def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
206                           dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 +  &
207                           dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
208
209                      IF ( def < 0.0 )  def = 0.0
210
211                      tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
212
213                   ENDDO
214
215                ENDIF
216
217             ENDDO
218
219!
220!--          (4) - will allways be executed.
221!--          'special case: free atmosphere' (as for case (0))
222             DO  j = nys, nyn
223
224                IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) &
225                THEN
226
227                   k = nzb_diff_s_outer(j,i)-1
228
229                   dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
230                   dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
231                                    u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
232                   dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
233                                    u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
234
235                   dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
236                                    v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
237                   dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
238                   dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
239                                    v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
240
241                   dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
242                                    w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
243                   dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
244                                    w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
245                   dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
246
247                   def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
248                         dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
249                         dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
250
251                   IF ( def < 0.0 )  def = 0.0
252
253                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
254
255                ENDIF
256
257             ENDDO
258
259!
260!--          Position without adjacent wall
261!--          (1) - will allways be executed.
262!--          'bottom only: use u_0,v_0'
263             DO  j = nys, nyn
264
265                IF ( ( wall_e_x(j,i) == 0.0 ) .AND. ( wall_e_y(j,i) == 0.0 ) ) &
266                THEN
267
268                   k = nzb_diff_s_inner(j,i)-1
269
270                   dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
271                   dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
272                                    u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
273                   dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
274                                    u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
275
276                   dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
277                                    v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
278                   dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
279                   dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
280                                    v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
281
282                   dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
283                                    w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
284                   dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
285                                    w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
286                   dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
287
288                   def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
289                         dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
290                         dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
291
292                   IF ( def < 0.0 )  def = 0.0
293
294                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
295               
296                ENDIF
297
298             ENDDO
299
300          ELSEIF ( use_surface_fluxes )  THEN
301
302             DO  j = nys, nyn
303
304                k = nzb_diff_s_outer(j,i)-1
305
306                dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
307                dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
308                                 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
309                dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
310                                 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
311
312                dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
313                                 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
314                dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
315                dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
316                                 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
317
318                dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
319                                 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
320                dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
321                                 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
322                dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
323
324                def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
325                      dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
326                      dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
327
328                IF ( def < 0.0 )  def = 0.0
329
330                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
331
332             ENDDO
333
334          ENDIF
335
336!
337!--       Calculate TKE production by buoyancy
338          IF ( .NOT. moisture )  THEN
339
340             DO  j = nys, nyn
341
342                DO  k = nzb_diff_s_inner(j,i), nzt_diff
343                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * &
344                                    ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
345                ENDDO
346
347                IF ( use_surface_fluxes )  THEN
348                   k = nzb_diff_s_inner(j,i)-1
349                   tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
350                ENDIF
351
352                IF ( use_top_fluxes )  THEN
353                   k = nzt
354                   tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
355                ENDIF
356
357             ENDDO
358
359          ELSE
360
361             DO  j = nys, nyn
362
363                DO  k = nzb_diff_s_inner(j,i), nzt_diff
364
365                   IF ( .NOT. cloud_physics )  THEN
366                      k1 = 1.0 + 0.61 * q(k,j,i)
367                      k2 = 0.61 * pt(k,j,i)
368                   ELSE
369                      IF ( ql(k,j,i) == 0.0 )  THEN
370                         k1 = 1.0 + 0.61 * q(k,j,i)
371                         k2 = 0.61 * pt(k,j,i)
372                      ELSE
373                         theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
374                         temp  = theta * t_d_pt(k)
375                         k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
376                                    ( q(k,j,i) - ql(k,j,i) ) *          &
377                              ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
378                              ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
379                              ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
380                         k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
381                      ENDIF
382                   ENDIF
383
384                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * &
385                                      ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
386                                        k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
387                                      ) * dd2zu(k)
388                ENDDO
389
390             ENDDO
391
392             IF ( use_surface_fluxes )  THEN
393
394                DO  j = nys, nyn
395
396                   k = nzb_diff_s_inner(j,i)-1
397
398                   IF ( .NOT. cloud_physics )  THEN
399                      k1 = 1.0 + 0.61 * q(k,j,i)
400                      k2 = 0.61 * pt(k,j,i)
401                   ELSE
402                      IF ( ql(k,j,i) == 0.0 )  THEN
403                         k1 = 1.0 + 0.61 * q(k,j,i)
404                         k2 = 0.61 * pt(k,j,i)
405                      ELSE
406                         theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
407                         temp  = theta * t_d_pt(k)
408                         k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
409                                    ( q(k,j,i) - ql(k,j,i) ) *          &
410                              ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
411                              ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
412                              ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
413                         k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
414                      ENDIF
415                   ENDIF
416
417                   tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
418                                         ( k1* shf(j,i) + k2 * qsws(j,i) )
419                ENDDO
420
421             ENDIF
422
423             IF ( use_top_fluxes )  THEN
424
425                DO  j = nys, nyn
426
427                   k = nzt
428
429                   IF ( .NOT. cloud_physics )  THEN
430                      k1 = 1.0 + 0.61 * q(k,j,i)
431                      k2 = 0.61 * pt(k,j,i)
432                   ELSE
433                      IF ( ql(k,j,i) == 0.0 )  THEN
434                         k1 = 1.0 + 0.61 * q(k,j,i)
435                         k2 = 0.61 * pt(k,j,i)
436                      ELSE
437                         theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
438                         temp  = theta * t_d_pt(k)
439                         k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
440                                    ( q(k,j,i) - ql(k,j,i) ) *          &
441                              ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
442                              ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
443                              ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
444                         k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
445                      ENDIF
446                   ENDIF
447
448                   tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
449                                         ( k1* tswst(j,i) + k2 * qswst(j,i) )
450                ENDDO
451
452             ENDIF
453
454          ENDIF
455
456       ENDDO
457
458    END SUBROUTINE production_e
459
460
461!------------------------------------------------------------------------------!
462! Call for grid point i,j
463!------------------------------------------------------------------------------!
464    SUBROUTINE production_e_ij( i, j )
465
466       USE arrays_3d
467       USE cloud_parameters
468       USE control_parameters
469       USE grid_variables
470       USE indices
471       USE statistics
472
473       IMPLICIT NONE
474
475       INTEGER ::  i, j, k
476
477       REAL    ::  def, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, &
478                   k1, k2, theta, temp
479
480       REAL, DIMENSION(nzb:nzt+1) ::  usvs, vsus, wsus,wsvs
481
482!
483!--    Calculate TKE production by shear
484       DO  k = nzb_diff_s_outer(j,i), nzt
485
486          dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
487          dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
488                           u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
489          dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
490                           u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
491
492          dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
493                           v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
494          dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
495          dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
496                           v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
497
498          dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
499                           w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
500          dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
501                           w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
502          dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
503
504          def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 )                       &
505                + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 &
506                + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
507
508          IF ( def < 0.0 )  def = 0.0
509
510          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
511               
512       ENDDO
513
514       IF ( prandtl_layer )  THEN
515
516          IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) )  THEN
517
518!
519!--          Position beneath wall
520!--          (2) - Will allways be executed.
521!--          'bottom and wall: use u_0,v_0 and wall functions'
522             k = nzb_diff_s_inner(j,i)-1
523
524             dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
525             dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
526                            u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
527             dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
528             dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
529                            v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
530             dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
531
532             IF ( wall_e_y(j,i) /= 0.0 )  THEN
533                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
534                                    usvs, 1.0, 0.0, 0.0, 0.0 )
535                dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
536                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
537                                    wsvs, 0.0, 0.0, 1.0, 0.0 )
538                dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
539             ELSE
540                dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
541                                u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
542                dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
543                                w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
544             ENDIF
545
546             IF ( wall_e_x(j,i) /= 0.0 )  THEN
547                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
548                                    vsus, 0.0, 1.0, 0.0, 0.0 )
549                dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
550                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
551                                    wsus, 0.0, 0.0, 0.0, 1.0 )
552                dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
553             ELSE
554                dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
555                                v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
556                dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
557                                w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
558             ENDIF
559
560             def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
561                   dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
562                   dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
563
564             IF ( def < 0.0 )  def = 0.0
565
566             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
567
568!
569!--          (3) - will be executed only, if there is at least one level
570!--          between (2) and (4), i.e. the topography must have a
571!--          minimum height of 2 dz. Wall fluxes for this case have
572!--          already been calculated for (2).
573!--          'wall only: use wall functions'
574             DO  k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2
575
576                dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
577                dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
578                               u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
579                dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
580                dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
581                               v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
582                dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
583
584                IF ( wall_e_y(j,i) /= 0.0 )  THEN
585                   dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
586                   dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
587                ELSE
588                   dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
589                                   u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
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
594                IF ( wall_e_x(j,i) /= 0.0 )  THEN
595                   dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
596                   dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
597                ELSE
598                   dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
599                                   v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
600                   dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
601                                   w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
602                ENDIF
603
604                def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
605                      dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
606                      dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
607
608                IF ( def < 0.0 )  def = 0.0
609
610                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
611
612             ENDDO
613
614!
615!--          (4) - will allways be executed.
616!--          'special case: free atmosphere' (as for case (0))
617             k = nzb_diff_s_outer(j,i)-1
618
619             dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
620             dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
621                              u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
622             dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
623                              u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
624
625             dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
626                              v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
627             dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
628             dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
629                              v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
630
631             dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
632                              w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
633             dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
634                              w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
635             dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
636
637             def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
638                   dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
639                   dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
640
641             IF ( def < 0.0 )  def = 0.0
642
643             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
644
645          ELSE
646
647!
648!--          Position without adjacent wall
649!--          (1) - will allways be executed.
650!--          'bottom only: use u_0,v_0'
651             k = nzb_diff_s_inner(j,i)-1
652
653             dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
654             dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
655                              u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
656             dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
657                              u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
658
659             dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
660                              v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
661             dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
662             dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
663                              v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
664
665             dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
666                              w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
667             dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
668                              w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
669             dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
670
671             def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 )                       &
672                   + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 &
673                   + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
674
675             IF ( def < 0.0 )  def = 0.0
676
677             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
678
679          ENDIF
680
681       ELSEIF ( use_surface_fluxes )  THEN
682
683          k = nzb_diff_s_outer(j,i)-1
684
685          dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
686          dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
687                           u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
688          dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
689                           u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
690
691          dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
692                           v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
693          dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
694          dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
695                           v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
696
697          dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
698                           w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
699          dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
700                           w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
701          dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
702
703          def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
704                dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
705                dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
706
707          IF ( def < 0.0 )  def = 0.0
708
709          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
710
711       ENDIF
712
713!
714!--    Calculate TKE production by buoyancy
715       IF ( .NOT. moisture )  THEN
716
717          DO  k = nzb_diff_s_inner(j,i), nzt_diff
718             tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * &
719                                    ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
720          ENDDO
721
722          IF ( use_surface_fluxes )  THEN
723             k = nzb_diff_s_inner(j,i)-1
724             tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
725          ENDIF
726
727          IF ( use_top_fluxes )  THEN
728             k = nzt
729             tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
730          ENDIF
731
732       ELSE
733
734          DO  k = nzb_diff_s_inner(j,i), nzt_diff
735
736             IF ( .NOT. cloud_physics )  THEN
737                k1 = 1.0 + 0.61 * q(k,j,i)
738                k2 = 0.61 * pt(k,j,i)
739             ELSE
740                IF ( ql(k,j,i) == 0.0 )  THEN
741                   k1 = 1.0 + 0.61 * q(k,j,i)
742                   k2 = 0.61 * pt(k,j,i)
743                ELSE
744                   theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
745                   temp  = theta * t_d_pt(k)
746                   k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
747                              ( q(k,j,i) - ql(k,j,i) ) *          &
748                        ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
749                        ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
750                        ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
751                   k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
752                ENDIF
753             ENDIF
754
755             tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *      &
756                                      ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
757                                        k2 * ( q(k+1,j,i) - q(k-1,j,i) )   &
758                                      ) * dd2zu(k)
759          ENDDO
760
761          IF ( use_surface_fluxes )  THEN
762             k = nzb_diff_s_inner(j,i)-1
763
764             IF ( .NOT. cloud_physics ) THEN
765                k1 = 1.0 + 0.61 * q(k,j,i)
766                k2 = 0.61 * pt(k,j,i)
767             ELSE
768                IF ( ql(k,j,i) == 0.0 )  THEN
769                   k1 = 1.0 + 0.61 * q(k,j,i)
770                   k2 = 0.61 * pt(k,j,i)
771                ELSE
772                   theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
773                   temp  = theta * t_d_pt(k)
774                   k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
775                              ( q(k,j,i) - ql(k,j,i) ) *          &
776                        ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
777                        ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
778                        ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
779                   k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
780                ENDIF
781             ENDIF
782
783             tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
784                                         ( k1* shf(j,i) + k2 * qsws(j,i) )
785          ENDIF
786
787          IF ( use_top_fluxes )  THEN
788             k = nzt
789
790             IF ( .NOT. cloud_physics ) THEN
791                k1 = 1.0 + 0.61 * q(k,j,i)
792                k2 = 0.61 * pt(k,j,i)
793             ELSE
794                IF ( ql(k,j,i) == 0.0 )  THEN
795                   k1 = 1.0 + 0.61 * q(k,j,i)
796                   k2 = 0.61 * pt(k,j,i)
797                ELSE
798                   theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
799                   temp  = theta * t_d_pt(k)
800                   k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
801                              ( q(k,j,i) - ql(k,j,i) ) *          &
802                        ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
803                        ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
804                        ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
805                   k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
806                ENDIF
807             ENDIF
808
809             tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
810                                         ( k1* tswst(j,i) + k2 * qswst(j,i) )
811          ENDIF
812
813       ENDIF
814
815    END SUBROUTINE production_e_ij
816
817
818    SUBROUTINE production_e_init
819
820       USE arrays_3d
821       USE control_parameters
822       USE grid_variables
823       USE indices
824
825       IMPLICIT NONE
826
827       INTEGER ::  i, j, ku, kv
828
829       IF ( prandtl_layer )  THEN
830
831          IF ( first_call )  THEN
832             ALLOCATE( u_0(nys-1:nyn+1,nxl-1:nxr+1), &
833                       v_0(nys-1:nyn+1,nxl-1:nxr+1) )
834             first_call = .FALSE.
835          ENDIF
836
837!
838!--       Calculate a virtual velocity at the surface in a way that the
839!--       vertical velocity gradient at k = 1 (u(k+1)-u_0) matches the
840!--       Prandtl law (-w'u'/km). This gradient is used in the TKE shear
841!--       production term at k=1 (see production_e_ij).
842!--       The velocity gradient has to be limited in case of too small km
843!--       (otherwise the timestep may be significantly reduced by large
844!--       surface winds).
845!--       WARNING: the exact analytical solution would require the determination
846!--                of the eddy diffusivity by km = u* * kappa * zp / phi_m.
847          !$OMP PARALLEL DO PRIVATE( ku, kv )
848          DO  i = nxl, nxr
849             DO  j = nys, nyn
850
851                ku = nzb_u_inner(j,i)+1
852                kv = nzb_v_inner(j,i)+1
853
854                u_0(j,i) = u(ku+1,j,i) + usws(j,i) * ( zu(ku+1) - zu(ku-1) ) / &
855                                 ( 0.5 * ( km(ku,j,i) + km(ku,j,i-1) ) +       &
856                                   1.0E-20 )
857!                                  ( us(j,i) * kappa * zu(1) )
858                v_0(j,i) = v(kv+1,j,i) + vsws(j,i) * ( zu(kv+1) - zu(kv-1) ) / &
859                                 ( 0.5 * ( km(kv,j,i) + km(kv,j-1,i) ) +       &
860                                   1.0E-20 )
861!                                  ( us(j,i) * kappa * zu(1) )
862
863                IF ( ABS( u(ku+1,j,i) - u_0(j,i) )  > &
864                     ABS( u(ku+1,j,i) - u(ku-1,j,i) ) )  u_0(j,i) = u(ku-1,j,i)
865                IF ( ABS( v(kv+1,j,i) - v_0(j,i) )  > &
866                     ABS( v(kv+1,j,i) - v(kv-1,j,i) ) )  v_0(j,i) = v(kv-1,j,i)
867
868             ENDDO
869          ENDDO
870
871          CALL exchange_horiz_2d( u_0 )
872          CALL exchange_horiz_2d( v_0 )
873
874       ENDIF
875
876    END SUBROUTINE production_e_init
877
878 END MODULE production_e_mod
Note: See TracBrowser for help on using the repository browser.