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

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

further checkin of preliminary changes

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