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

Last change on this file since 187 was 187, checked in by letzel, 13 years ago
  • new: descriptions of plant canopy model and prandtl layer (trunk/DOC/misc)
  • changed: more consistent flux definitions; modification of the integrated version of the profile function for momentum for unstable stratification (wall_fluxes, production_e)
  • bugfix: change definition of us_wall from 1D to 2D (prandtl_fluxes, wall_fluxes)
  • Property svn:keywords set to Id
File size: 40.6 KB
Line 
1 MODULE production_e_mod
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6! Change: add 'minus' sign to fluxes obtained from subroutine wall_fluxes_e for
7!         consistency with subroutine wall_fluxes
8!
9! Former revisions:
10! -----------------
11! $Id: production_e.f90 187 2008-08-06 16:25:09Z letzel $
12!
13! 124 2007-10-19 15:47:46Z raasch
14! Bugfix: calculation of density flux in the ocean now starts from nzb+1
15!
16! 108 2007-08-24 15:10:38Z letzel
17! Bugfix: wrong sign removed from the buoyancy production term in the case
18! use_reference = .T.,
19! u_0 and v_0 are calculated for nxr+1, nyn+1 also (otherwise these values are
20! not available in case of non-cyclic boundary conditions)
21! Bugfix for ocean density flux at bottom
22!
23! 97 2007-06-21 08:23:15Z raasch
24! Energy production by density flux (in ocean) added
25! use_pt_reference renamed use_reference
26!
27! 75 2007-03-22 09:54:05Z raasch
28! Wall functions now include diabatic conditions, call of routine wall_fluxes_e,
29! reference temperature pt_reference can be used in buoyancy term,
30! moisture renamed humidity
31!
32! 37 2007-03-01 08:33:54Z raasch
33! Calculation extended for gridpoint nzt, extended for given temperature /
34! humidity fluxes at the top, wall-part is now executed in case that a
35! Prandtl-layer is switched on (instead of surfaces fluxes switched on)
36!
37! RCS Log replace by Id keyword, revision history cleaned up
38!
39! Revision 1.21  2006/04/26 12:45:35  raasch
40! OpenMP parallelization of production_e_init
41!
42! Revision 1.1  1997/09/19 07:45:35  raasch
43! Initial revision
44!
45!
46! Description:
47! ------------
48! Production terms (shear + buoyancy) of the TKE
49! WARNING: the case with prandtl_layer = F and use_surface_fluxes = T is
50!          not considered well!
51!------------------------------------------------------------------------------!
52
53    USE wall_fluxes_mod
54
55    PRIVATE
56    PUBLIC production_e, production_e_init
57
58    LOGICAL, SAVE ::  first_call = .TRUE.
59
60    REAL, DIMENSION(:,:), ALLOCATABLE, SAVE ::  u_0, v_0
61
62    INTERFACE production_e
63       MODULE PROCEDURE production_e
64       MODULE PROCEDURE production_e_ij
65    END INTERFACE production_e
66   
67    INTERFACE production_e_init
68       MODULE PROCEDURE production_e_init
69    END INTERFACE production_e_init
70 
71 CONTAINS
72
73
74!------------------------------------------------------------------------------!
75! Call for all grid points
76!------------------------------------------------------------------------------!
77    SUBROUTINE production_e
78
79       USE arrays_3d
80       USE cloud_parameters
81       USE control_parameters
82       USE grid_variables
83       USE indices
84       USE statistics
85
86       IMPLICIT NONE
87
88       INTEGER ::  i, j, k
89
90       REAL    ::  def, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, &
91                   k1, k2, theta, temp
92
93!       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  usvs, vsus, wsus, wsvs
94       REAL, DIMENSION(nzb:nzt+1) ::   usvs, vsus, wsus, wsvs
95
96!
97!--    First calculate horizontal momentum flux u'v', w'v', v'u', w'u' at
98!--    vertical walls, if neccessary
99!--    So far, results are slightly different from the ij-Version.
100!--    Therefore, ij-Version is called further below within the ij-loops.
101!       IF ( topography /= 'flat' )  THEN
102!          CALL wall_fluxes_e( usvs, 1.0, 0.0, 0.0, 0.0, wall_e_y )
103!          CALL wall_fluxes_e( wsvs, 0.0, 0.0, 1.0, 0.0, wall_e_y )
104!          CALL wall_fluxes_e( vsus, 0.0, 1.0, 0.0, 0.0, wall_e_x )
105!          CALL wall_fluxes_e( wsus, 0.0, 0.0, 0.0, 1.0, wall_e_x )
106!       ENDIF
107
108!
109!--    Calculate TKE production by shear
110       DO  i = nxl, nxr
111
112          DO  j = nys, nyn
113             DO  k = nzb_diff_s_outer(j,i), nzt
114
115                dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
116                dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
117                                 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
118                dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
119                                 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
120
121                dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
122                                 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
123                dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
124                dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
125                                 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
126
127                dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
128                                 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
129                dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
130                                 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
131                dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
132
133                def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
134                      dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
135                      dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
136
137                IF ( def < 0.0 )  def = 0.0
138
139                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
140               
141             ENDDO
142          ENDDO
143
144          IF ( prandtl_layer )  THEN
145
146!
147!--          Position beneath wall
148!--          (2) - Will allways be executed.
149!--          'bottom and wall: use u_0,v_0 and wall functions'
150             DO  j = nys, nyn
151
152                IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) &
153                THEN
154
155                   k = nzb_diff_s_inner(j,i) - 1
156                   dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
157                   dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
158                                  u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
159                   dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
160                   dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
161                                  v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
162                   dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
163
164                   IF ( wall_e_y(j,i) /= 0.0 )  THEN
165                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
166                                          usvs, 1.0, 0.0, 0.0, 0.0 )
167                      dudy = - wall_e_y(j,i) * usvs(k) / km(k,j,i)
168!                      dudy = - wall_e_y(j,i) * usvs(k,j,i) / km(k,j,i)
169                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
170                                          wsvs, 0.0, 0.0, 1.0, 0.0 )
171                      dwdy = - wall_e_y(j,i) * wsvs(k) / km(k,j,i)
172!                      dwdy = - wall_e_y(j,i) * wsvs(k,j,i) / km(k,j,i)
173                   ELSE
174                      dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
175                                      u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
176                      dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
177                                      w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
178                   ENDIF
179
180                   IF ( wall_e_x(j,i) /= 0.0 )  THEN
181                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
182                                          vsus, 0.0, 1.0, 0.0, 0.0 )
183                      dvdx = - wall_e_x(j,i) * vsus(k) / km(k,j,i)
184!                      dvdx = - wall_e_x(j,i) * vsus(k,j,i) / km(k,j,i)
185                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
186                                          wsus, 0.0, 0.0, 0.0, 1.0 )
187                      dwdx = - wall_e_x(j,i) * wsus(k) / km(k,j,i)
188!                      dwdx = - wall_e_x(j,i) * wsus(k,j,i) / km(k,j,i)
189                   ELSE
190                      dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
191                                      v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
192                      dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
193                                      w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
194                   ENDIF
195
196                   def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
197                         dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
198                         dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
199
200                   IF ( def < 0.0 )  def = 0.0
201
202                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
203
204
205!
206!--                (3) - will be executed only, if there is at least one level
207!--                between (2) and (4), i.e. the topography must have a
208!--                minimum height of 2 dz. Wall fluxes for this case have
209!--                already been calculated for (2).
210!--                'wall only: use wall functions'
211
212                   DO  k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2
213
214                      dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
215                      dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
216                                     u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
217                      dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
218                      dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
219                                     v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
220                      dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
221
222                      IF ( wall_e_y(j,i) /= 0.0 )  THEN
223                         dudy = - wall_e_y(j,i) * usvs(k) / km(k,j,i)
224!                         dudy = - wall_e_y(j,i) * usvs(k,j,i) / km(k,j,i)
225                         dwdy = - wall_e_y(j,i) * wsvs(k) / km(k,j,i)
226!                         dwdy = - wall_e_y(j,i) * wsvs(k,j,i) / km(k,j,i)
227                      ELSE
228                         dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
229                                         u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
230                         dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
231                                         w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
232                      ENDIF
233
234                      IF ( wall_e_x(j,i) /= 0.0 )  THEN
235                         dvdx = - wall_e_x(j,i) * vsus(k) / km(k,j,i)
236!                         dvdx = - wall_e_x(j,i) * vsus(k,j,i) / km(k,j,i)
237                         dwdx = - wall_e_x(j,i) * wsus(k) / km(k,j,i)
238!                         dwdx = - wall_e_x(j,i) * wsus(k,j,i) / km(k,j,i)
239                      ELSE
240                         dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
241                                         v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
242                         dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
243                                         w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
244                      ENDIF
245
246                      def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
247                           dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 +  &
248                           dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
249
250                      IF ( def < 0.0 )  def = 0.0
251
252                      tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
253
254                   ENDDO
255
256                ENDIF
257
258             ENDDO
259
260!
261!--          (4) - will allways be executed.
262!--          'special case: free atmosphere' (as for case (0))
263             DO  j = nys, nyn
264
265                IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) &
266                THEN
267
268                   k = nzb_diff_s_outer(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(k-1,j,i) - u(k-1,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(k-1,j,i) - v(k-1,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!
301!--          Position without adjacent wall
302!--          (1) - will allways be executed.
303!--          'bottom only: use u_0,v_0'
304             DO  j = nys, nyn
305
306                IF ( ( wall_e_x(j,i) == 0.0 ) .AND. ( wall_e_y(j,i) == 0.0 ) ) &
307                THEN
308
309                   k = nzb_diff_s_inner(j,i)-1
310
311                   dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
312                   dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
313                                    u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
314                   dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
315                                    u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
316
317                   dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
318                                    v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
319                   dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
320                   dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
321                                    v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
322
323                   dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
324                                    w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
325                   dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
326                                    w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
327                   dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
328
329                   def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
330                         dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
331                         dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
332
333                   IF ( def < 0.0 )  def = 0.0
334
335                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
336               
337                ENDIF
338
339             ENDDO
340
341          ELSEIF ( use_surface_fluxes )  THEN
342
343             DO  j = nys, nyn
344
345                k = nzb_diff_s_outer(j,i)-1
346
347                dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
348                dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
349                                 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
350                dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
351                                 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
352
353                dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
354                                 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
355                dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
356                dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
357                                 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
358
359                dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
360                                 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
361                dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
362                                 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
363                dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
364
365                def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
366                      dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
367                      dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
368
369                IF ( def < 0.0 )  def = 0.0
370
371                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
372
373             ENDDO
374
375          ENDIF
376
377!
378!--       Calculate TKE production by buoyancy
379          IF ( .NOT. humidity )  THEN
380
381             IF ( use_reference )  THEN
382
383                IF ( ocean )  THEN
384!
385!--                So far in the ocean no special treatment of density flux in
386!--                the bottom and top surface layer
387                   DO  j = nys, nyn
388                      DO  k = nzb_s_inner(j,i)+1, nzt
389                         tend(k,j,i) = tend(k,j,i) +                    &
390                                       kh(k,j,i) * g / prho_reference * &
391                                       ( rho(k+1,j,i)-rho(k-1,j,i) ) * dd2zu(k)
392                      ENDDO
393                   ENDDO
394
395                ELSE
396
397                   DO  j = nys, nyn
398                      DO  k = nzb_diff_s_inner(j,i), nzt_diff
399                         tend(k,j,i) = tend(k,j,i) -                  &
400                                       kh(k,j,i) * g / pt_reference * &
401                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
402                      ENDDO
403
404                      IF ( use_surface_fluxes )  THEN
405                         k = nzb_diff_s_inner(j,i)-1
406                         tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
407                      ENDIF
408
409                      IF ( use_top_fluxes )  THEN
410                         k = nzt
411                         tend(k,j,i) = tend(k,j,i) + g / pt_reference * &
412                                                     tswst(j,i)
413                      ENDIF
414                   ENDDO
415
416                ENDIF
417
418             ELSE
419
420                IF ( ocean )  THEN
421!
422!--                So far in the ocean no special treatment of density flux in
423!--                the bottom and top surface layer
424                   DO  j = nys, nyn
425                      DO  k = nzb_s_inner(j,i)+1, nzt
426                         tend(k,j,i) = tend(k,j,i) -                &
427                                       kh(k,j,i) * g / rho(k,j,i) * &
428                                       ( rho(k+1,j,i)-rho(k-1,j,i) ) * dd2zu(k)
429                      ENDDO
430                   ENDDO
431
432                ELSE
433
434                   DO  j = nys, nyn
435                      DO  k = nzb_diff_s_inner(j,i), nzt_diff
436                         tend(k,j,i) = tend(k,j,i) -               &
437                                       kh(k,j,i) * g / pt(k,j,i) * &
438                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
439                      ENDDO
440
441                      IF ( use_surface_fluxes )  THEN
442                         k = nzb_diff_s_inner(j,i)-1
443                         tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
444                      ENDIF
445
446                      IF ( use_top_fluxes )  THEN
447                         k = nzt
448                         tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
449                      ENDIF
450                   ENDDO
451
452                ENDIF
453
454             ENDIF
455
456          ELSE
457
458             DO  j = nys, nyn
459
460                DO  k = nzb_diff_s_inner(j,i), nzt_diff
461
462                   IF ( .NOT. cloud_physics )  THEN
463                      k1 = 1.0 + 0.61 * q(k,j,i)
464                      k2 = 0.61 * pt(k,j,i)
465                   ELSE
466                      IF ( ql(k,j,i) == 0.0 )  THEN
467                         k1 = 1.0 + 0.61 * q(k,j,i)
468                         k2 = 0.61 * pt(k,j,i)
469                      ELSE
470                         theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
471                         temp  = theta * t_d_pt(k)
472                         k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
473                                    ( q(k,j,i) - ql(k,j,i) ) *          &
474                              ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
475                              ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
476                              ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
477                         k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
478                      ENDIF
479                   ENDIF
480
481                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * &
482                                      ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
483                                        k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
484                                      ) * dd2zu(k)
485                ENDDO
486
487             ENDDO
488
489             IF ( use_surface_fluxes )  THEN
490
491                DO  j = nys, nyn
492
493                   k = nzb_diff_s_inner(j,i)-1
494
495                   IF ( .NOT. cloud_physics )  THEN
496                      k1 = 1.0 + 0.61 * q(k,j,i)
497                      k2 = 0.61 * pt(k,j,i)
498                   ELSE
499                      IF ( ql(k,j,i) == 0.0 )  THEN
500                         k1 = 1.0 + 0.61 * q(k,j,i)
501                         k2 = 0.61 * pt(k,j,i)
502                      ELSE
503                         theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
504                         temp  = theta * t_d_pt(k)
505                         k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
506                                    ( q(k,j,i) - ql(k,j,i) ) *          &
507                              ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
508                              ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
509                              ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
510                         k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
511                      ENDIF
512                   ENDIF
513
514                   tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
515                                         ( k1* shf(j,i) + k2 * qsws(j,i) )
516                ENDDO
517
518             ENDIF
519
520             IF ( use_top_fluxes )  THEN
521
522                DO  j = nys, nyn
523
524                   k = nzt
525
526                   IF ( .NOT. cloud_physics )  THEN
527                      k1 = 1.0 + 0.61 * q(k,j,i)
528                      k2 = 0.61 * pt(k,j,i)
529                   ELSE
530                      IF ( ql(k,j,i) == 0.0 )  THEN
531                         k1 = 1.0 + 0.61 * q(k,j,i)
532                         k2 = 0.61 * pt(k,j,i)
533                      ELSE
534                         theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
535                         temp  = theta * t_d_pt(k)
536                         k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
537                                    ( q(k,j,i) - ql(k,j,i) ) *          &
538                              ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
539                              ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
540                              ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
541                         k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
542                      ENDIF
543                   ENDIF
544
545                   tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
546                                         ( k1* tswst(j,i) + k2 * qswst(j,i) )
547                ENDDO
548
549             ENDIF
550
551          ENDIF
552
553       ENDDO
554
555    END SUBROUTINE production_e
556
557
558!------------------------------------------------------------------------------!
559! Call for grid point i,j
560!------------------------------------------------------------------------------!
561    SUBROUTINE production_e_ij( i, j )
562
563       USE arrays_3d
564       USE cloud_parameters
565       USE control_parameters
566       USE grid_variables
567       USE indices
568       USE statistics
569
570       IMPLICIT NONE
571
572       INTEGER ::  i, j, k
573
574       REAL    ::  def, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, &
575                   k1, k2, theta, temp
576
577       REAL, DIMENSION(nzb:nzt+1) ::  usvs, vsus, wsus, wsvs
578
579!
580!--    Calculate TKE production by shear
581       DO  k = nzb_diff_s_outer(j,i), nzt
582
583          dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
584          dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
585                           u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
586          dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
587                           u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
588
589          dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
590                           v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
591          dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
592          dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
593                           v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
594
595          dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
596                           w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
597          dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
598                           w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
599          dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
600
601          def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 )                       &
602                + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 &
603                + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
604
605          IF ( def < 0.0 )  def = 0.0
606
607          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
608               
609       ENDDO
610
611       IF ( prandtl_layer )  THEN
612
613          IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) )  THEN
614
615!
616!--          Position beneath wall
617!--          (2) - Will allways be executed.
618!--          'bottom and wall: use u_0,v_0 and wall functions'
619             k = nzb_diff_s_inner(j,i)-1
620
621             dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
622             dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
623                            u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
624             dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
625             dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
626                            v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
627             dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
628
629             IF ( wall_e_y(j,i) /= 0.0 )  THEN
630                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
631                                    usvs, 1.0, 0.0, 0.0, 0.0 )
632                dudy = - wall_e_y(j,i) * usvs(k) / km(k,j,i)
633                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
634                                    wsvs, 0.0, 0.0, 1.0, 0.0 )
635                dwdy = - wall_e_y(j,i) * wsvs(k) / km(k,j,i)
636             ELSE
637                dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
638                                u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
639                dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
640                                w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
641             ENDIF
642
643             IF ( wall_e_x(j,i) /= 0.0 )  THEN
644                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
645                                    vsus, 0.0, 1.0, 0.0, 0.0 )
646                dvdx = - wall_e_x(j,i) * vsus(k) / km(k,j,i)
647                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
648                                    wsus, 0.0, 0.0, 0.0, 1.0 )
649                dwdx = - wall_e_x(j,i) * wsus(k) / km(k,j,i)
650             ELSE
651                dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
652                                v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
653                dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
654                                w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
655             ENDIF
656
657             def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
658                   dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
659                   dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
660
661             IF ( def < 0.0 )  def = 0.0
662
663             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
664
665!
666!--          (3) - will be executed only, if there is at least one level
667!--          between (2) and (4), i.e. the topography must have a
668!--          minimum height of 2 dz. Wall fluxes for this case have
669!--          already been calculated for (2).
670!--          'wall only: use wall functions'
671             DO  k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2
672
673                dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
674                dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
675                               u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
676                dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
677                dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
678                               v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
679                dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
680
681                IF ( wall_e_y(j,i) /= 0.0 )  THEN
682                   dudy = - wall_e_y(j,i) * usvs(k) / km(k,j,i)
683                   dwdy = - wall_e_y(j,i) * wsvs(k) / km(k,j,i)
684                ELSE
685                   dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
686                                   u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
687                   dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
688                                   w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
689                ENDIF
690
691                IF ( wall_e_x(j,i) /= 0.0 )  THEN
692                   dvdx = - wall_e_x(j,i) * vsus(k) / km(k,j,i)
693                   dwdx = - wall_e_x(j,i) * wsus(k) / km(k,j,i)
694                ELSE
695                   dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
696                                   v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
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                ENDIF
700
701                def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
702                      dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
703                      dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
704
705                IF ( def < 0.0 )  def = 0.0
706
707                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
708
709             ENDDO
710
711!
712!--          (4) - will allways be executed.
713!--          'special case: free atmosphere' (as for case (0))
714             k = nzb_diff_s_outer(j,i)-1
715
716             dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
717             dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
718                              u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
719             dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
720                              u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
721
722             dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
723                              v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
724             dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
725             dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
726                              v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
727
728             dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
729                              w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
730             dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
731                              w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
732             dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
733
734             def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
735                   dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
736                   dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
737
738             IF ( def < 0.0 )  def = 0.0
739
740             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
741
742          ELSE
743
744!
745!--          Position without adjacent wall
746!--          (1) - will allways be executed.
747!--          'bottom only: use u_0,v_0'
748             k = nzb_diff_s_inner(j,i)-1
749
750             dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
751             dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
752                              u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
753             dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
754                              u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
755
756             dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
757                              v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
758             dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
759             dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
760                              v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
761
762             dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
763                              w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
764             dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
765                              w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
766             dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
767
768             def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 )                       &
769                   + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 &
770                   + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
771
772             IF ( def < 0.0 )  def = 0.0
773
774             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
775
776          ENDIF
777
778       ELSEIF ( use_surface_fluxes )  THEN
779
780          k = nzb_diff_s_outer(j,i)-1
781
782          dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
783          dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
784                           u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
785          dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
786                           u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
787
788          dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
789                           v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
790          dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
791          dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
792                           v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
793
794          dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
795                           w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
796          dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
797                           w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
798          dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
799
800          def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
801                dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
802                dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
803
804          IF ( def < 0.0 )  def = 0.0
805
806          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
807
808       ENDIF
809
810!
811!--    Calculate TKE production by buoyancy
812       IF ( .NOT. humidity )  THEN
813
814          IF ( use_reference )  THEN
815
816             IF ( ocean )  THEN
817!
818!--             So far in the ocean no special treatment of density flux in the
819!--             bottom and top surface layer
820                DO  k = nzb_s_inner(j,i)+1, nzt
821                   tend(k,j,i) = tend(k,j,i) + kh(k,j,i) * g / prho_reference * &
822                                      ( rho(k+1,j,i) - rho(k-1,j,i) ) * dd2zu(k)
823                ENDDO
824
825             ELSE
826
827                DO  k = nzb_diff_s_inner(j,i), nzt_diff
828                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt_reference * &
829                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
830                ENDDO
831
832                IF ( use_surface_fluxes )  THEN
833                   k = nzb_diff_s_inner(j,i)-1
834                   tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
835                ENDIF
836
837                IF ( use_top_fluxes )  THEN
838                   k = nzt
839                   tend(k,j,i) = tend(k,j,i) + g / pt_reference * tswst(j,i)
840                ENDIF
841
842             ENDIF
843
844          ELSE
845
846             IF ( ocean )  THEN
847!
848!--             So far in the ocean no special treatment of density flux in the
849!--             bottom and top surface layer
850                DO  k = nzb_s_inner(j,i)+1, nzt
851                   tend(k,j,i) = tend(k,j,i) + kh(k,j,i) * g / rho(k,j,i) * &
852                                      ( rho(k+1,j,i) - rho(k-1,j,i) ) * dd2zu(k)
853                ENDDO
854
855             ELSE
856
857                DO  k = nzb_diff_s_inner(j,i), nzt_diff
858                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * &
859                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
860                ENDDO
861
862                IF ( use_surface_fluxes )  THEN
863                   k = nzb_diff_s_inner(j,i)-1
864                   tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
865                ENDIF
866
867                IF ( use_top_fluxes )  THEN
868                   k = nzt
869                   tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
870                ENDIF
871
872             ENDIF
873
874          ENDIF
875
876       ELSE
877
878          DO  k = nzb_diff_s_inner(j,i), nzt_diff
879
880             IF ( .NOT. cloud_physics )  THEN
881                k1 = 1.0 + 0.61 * q(k,j,i)
882                k2 = 0.61 * pt(k,j,i)
883             ELSE
884                IF ( ql(k,j,i) == 0.0 )  THEN
885                   k1 = 1.0 + 0.61 * q(k,j,i)
886                   k2 = 0.61 * pt(k,j,i)
887                ELSE
888                   theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
889                   temp  = theta * t_d_pt(k)
890                   k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
891                              ( q(k,j,i) - ql(k,j,i) ) *          &
892                        ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
893                        ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
894                        ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
895                   k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
896                ENDIF
897             ENDIF
898
899             tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *      &
900                                      ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
901                                        k2 * ( q(k+1,j,i) - q(k-1,j,i) )   &
902                                      ) * dd2zu(k)
903          ENDDO
904
905          IF ( use_surface_fluxes )  THEN
906             k = nzb_diff_s_inner(j,i)-1
907
908             IF ( .NOT. cloud_physics ) THEN
909                k1 = 1.0 + 0.61 * q(k,j,i)
910                k2 = 0.61 * pt(k,j,i)
911             ELSE
912                IF ( ql(k,j,i) == 0.0 )  THEN
913                   k1 = 1.0 + 0.61 * q(k,j,i)
914                   k2 = 0.61 * pt(k,j,i)
915                ELSE
916                   theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
917                   temp  = theta * t_d_pt(k)
918                   k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
919                              ( q(k,j,i) - ql(k,j,i) ) *          &
920                        ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
921                        ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
922                        ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
923                   k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
924                ENDIF
925             ENDIF
926
927             tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
928                                         ( k1* shf(j,i) + k2 * qsws(j,i) )
929          ENDIF
930
931          IF ( use_top_fluxes )  THEN
932             k = nzt
933
934             IF ( .NOT. cloud_physics ) THEN
935                k1 = 1.0 + 0.61 * q(k,j,i)
936                k2 = 0.61 * pt(k,j,i)
937             ELSE
938                IF ( ql(k,j,i) == 0.0 )  THEN
939                   k1 = 1.0 + 0.61 * q(k,j,i)
940                   k2 = 0.61 * pt(k,j,i)
941                ELSE
942                   theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
943                   temp  = theta * t_d_pt(k)
944                   k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
945                              ( q(k,j,i) - ql(k,j,i) ) *          &
946                        ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
947                        ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
948                        ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
949                   k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
950                ENDIF
951             ENDIF
952
953             tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
954                                         ( k1* tswst(j,i) + k2 * qswst(j,i) )
955          ENDIF
956
957       ENDIF
958
959    END SUBROUTINE production_e_ij
960
961
962    SUBROUTINE production_e_init
963
964       USE arrays_3d
965       USE control_parameters
966       USE grid_variables
967       USE indices
968
969       IMPLICIT NONE
970
971       INTEGER ::  i, j, ku, kv
972
973       IF ( prandtl_layer )  THEN
974
975          IF ( first_call )  THEN
976             ALLOCATE( u_0(nys-1:nyn+1,nxl-1:nxr+1), &
977                       v_0(nys-1:nyn+1,nxl-1:nxr+1) )
978             first_call = .FALSE.
979          ENDIF
980
981!
982!--       Calculate a virtual velocity at the surface in a way that the
983!--       vertical velocity gradient at k = 1 (u(k+1)-u_0) matches the
984!--       Prandtl law (-w'u'/km). This gradient is used in the TKE shear
985!--       production term at k=1 (see production_e_ij).
986!--       The velocity gradient has to be limited in case of too small km
987!--       (otherwise the timestep may be significantly reduced by large
988!--       surface winds).
989!--       Upper bounds are nxr+1 and nyn+1 because otherwise these values are
990!--       not available in case of non-cyclic boundary conditions.
991!--       WARNING: the exact analytical solution would require the determination
992!--                of the eddy diffusivity by km = u* * kappa * zp / phi_m.
993          !$OMP PARALLEL DO PRIVATE( ku, kv )
994          DO  i = nxl, nxr+1
995             DO  j = nys, nyn+1
996
997                ku = nzb_u_inner(j,i)+1
998                kv = nzb_v_inner(j,i)+1
999
1000                u_0(j,i) = u(ku+1,j,i) + usws(j,i) * ( zu(ku+1) - zu(ku-1) ) / &
1001                                 ( 0.5 * ( km(ku,j,i) + km(ku,j,i-1) ) +       &
1002                                   1.0E-20 )
1003!                                  ( us(j,i) * kappa * zu(1) )
1004                v_0(j,i) = v(kv+1,j,i) + vsws(j,i) * ( zu(kv+1) - zu(kv-1) ) / &
1005                                 ( 0.5 * ( km(kv,j,i) + km(kv,j-1,i) ) +       &
1006                                   1.0E-20 )
1007!                                  ( us(j,i) * kappa * zu(1) )
1008
1009                IF ( ABS( u(ku+1,j,i) - u_0(j,i) )  > &
1010                     ABS( u(ku+1,j,i) - u(ku-1,j,i) ) )  u_0(j,i) = u(ku-1,j,i)
1011                IF ( ABS( v(kv+1,j,i) - v_0(j,i) )  > &
1012                     ABS( v(kv+1,j,i) - v(kv-1,j,i) ) )  v_0(j,i) = v(kv-1,j,i)
1013
1014             ENDDO
1015          ENDDO
1016
1017          CALL exchange_horiz_2d( u_0 )
1018          CALL exchange_horiz_2d( v_0 )
1019
1020       ENDIF
1021
1022    END SUBROUTINE production_e_init
1023
1024 END MODULE production_e_mod
Note: See TracBrowser for help on using the repository browser.