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

Last change on this file since 759 was 759, checked in by raasch, 11 years ago

New:
---

The number of parallel I/O operations can be limited with new mrun-option -w.
(advec_particles, data_output_2d, data_output_3d, header, init_grid, init_pegrid, init_3d_model, modules, palm, parin, write_3d_binary)

Changed:


mrun option -T is obligatory

Errors:


Bugfix: No zero assignments to volume_flow_initial and volume_flow_area in
case of normal restart runs. (init_3d_model)

initialization of u_0, v_0. This is just to avoid access of uninitialized
memory in exchange_horiz_2d, which causes respective error messages
when the Intel thread checker (inspector) is used. (production_e)

Bugfix for ts limitation (prandtl_fluxes)

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