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

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

updating comments and rc-file

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