Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

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

Last change on this file since 97 was 97, checked in by raasch, 17 years ago

New:
---
ocean version including prognostic equation for salinity and equation of state for seawater. Routine buoyancy can be used with both temperature and density.
+ inipar-parameters bc_sa_t, bottom_salinityflux, ocean, sa_surface, sa_vertical_gradient, sa_vertical_gradient_level, top_salinityflux

advec_s_bc, average_3d_data, boundary_conds, buoyancy, check_parameters, data_output_2d, data_output_3d, diffusion_e, flow_statistics, header, init_grid, init_3d_model, modules, netcdf, parin, production_e, prognostic_equations, read_var_list, sum_up_3d_data, swap_timelevel, time_integration, user_interface, write_var_list, write_3d_binary

New:
eqn_state_seawater, init_ocean

Changed:

inipar-parameter use_pt_reference renamed use_reference

hydro_press renamed hyp, routine calc_mean_pt_profile renamed calc_mean_profile

format adjustments for the ocean version (run_control)

advec_particles, buoyancy, calc_liquid_water_content, check_parameters, diffusion_e, diffusivities, header, init_cloud_physics, modules, production_e, prognostic_equations, run_control

Errors:

Bugfix: height above topography instead of height above level k=0 is used for calculating the mixing length (diffusion_e and diffusivities).

Bugfix: error in boundary condition for TKE removed (advec_s_bc)

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