Home


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

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

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

New:
---

Plant canopy model of Watanabe (2004,BLM 112,307-341) added.
It can be switched on by the inipar parameter plant_canopy.
The inipar parameter canopy_mode can be used to prescribe a
plant canopy type. The default case is a homogeneous plant
canopy. Heterogeneous distributions of the leaf area
density and the canopy drag coefficient can be defined in the
new routine user_init_plant_canopy (user_interface).
The inipar parameters lad_surface, lad_vertical_gradient and
lad_vertical_gradient_level can be used in order to
prescribe the vertical profile of leaf area density. The
inipar parameter drag_coefficient determines the canopy
drag coefficient.
Finally, the inipar parameter pch_index determines the
index of the upper boundary of the plant canopy.

Allow new case bc_uv_t = 'dirichlet_0' for channel flow.

For unknown variables (CASE DEFAULT) call new subroutine user_data_output_dvrp

Pressure boundary conditions for vertical walls added to the multigrid solver.
They are applied using new wall flag arrays (wall_flags_..) which are defined
for each grid level. New argument gls added to routine user_init_grid
(user_interface).

Frequence of sorting particles can be controlled with new particles_par
parameter dt_sort_particles. Sorting is moved from the SGS timestep loop in
advec_particles after the end of this loop.

advec_particles, check_parameters, data_output_dvrp, header, init_3d_model, init_grid, init_particles, init_pegrid, modules, package_parin, parin, plant_canopy_model, read_var_list, read_3d_binary, user_interface, write_var_list, write_3d_binary

Changed:

Redefine initial nzb_local as the actual total size of topography (later the
extent of topography in nzb_local is reduced by 1dx at the E topography walls
and by 1dy at the N topography walls to form the basis for nzb_s_inner);
for consistency redefine 'single_building' case.

Vertical profiles now based on nzb_s_inner; they are divided by
ngp_2dh_s_inner (scalars, procucts of scalars) and ngp_2dh (staggered velocity
components and their products, procucts of scalars and velocity components),
respectively.

Allow two instead of one digit to specify isosurface and slicer variables.

Status of 3D-volume NetCDF data file only depends on switch netcdf_64bit_3d (check_open)

prognostic_equations include the respective wall_*flux in the parameter list of
calls of diffusion_s. Same as before, only the values of wall_heatflux(0:4)
can be assigned. At present, wall_humidityflux, wall_qflux, wall_salinityflux,
and wall_scalarflux are kept zero. diffusion_s uses the respective wall_*flux
instead of wall_heatflux. This update serves two purposes:

  • it avoids errors in calculations with humidity/scalar/salinity and prescribed

non-zero wall_heatflux,

  • it prepares PALM for a possible assignment of wall fluxes of

humidity/scalar/salinity in a future release.

buoyancy, check_open, data_output_dvrp, diffusion_s, diffusivities, flow_statistics, header, init_3d_model, init_dvrp, init_grid, modules, prognostic_equations

Errors:

Bugfix: summation of sums_l_l in diffusivities.

Several bugfixes in the ocean part: Initial density rho is calculated
(init_ocean). Error in initializing u_init and v_init removed
(check_parameters). Calculation of density flux now starts from
nzb+1 (production_e).

Bugfix: pleft/pright changed to pnorth/psouth in sendrecv of particle tail
numbers along y, small bugfixes in the SGS part (advec_particles)

Bugfix: model_string needed a default value (combine_plot_fields)

Bugfix: wavenumber calculation for even nx in routines maketri (poisfft)

Bugfix: assignment of fluxes at walls

Bugfix: absolute value of f must be used when calculating the Blackadar mixing length (init_1d_model)

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