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

Last change on this file since 106 was 106, checked in by raasch, 14 years ago

preliminary update of bugfixes and extensions for non-cyclic BCs

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