Home


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

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

Last change on this file since 402 was 392, checked in by raasch, 15 years ago

New:
---

Adapted for machine lck
(mrun, mbuild, subjob)

bc_lr/bc_ns in most subroutines replaced by LOGICAL variables bc_lr_cyc,
bc_ns_cyc for speed optimization
(check_parameters, diffusion_u, diffusion_v, diffusion_w, modules)

Additional timestep criterion in case of simulations with plant canopy (timestep)

Check for illegal entries in section_xy|xz|yz that exceed nz+1|ny+1|nx+1
(check_parameters)

Clipping of dvrp output implemented. Default colourtable for particles
implemented, particle attributes (color, dvrp_size) can be set with new
parameters particle_color, particle_dvrpsize, color_interval,
dvrpsize_interval (init_dvrp, data_output_dvrp, modules, user_data_output_dvrp).
Slicer attributes (dvrp) are set with new routine set_slicer_attributes_dvrp
and are controlled with existing parameters slicer_range_limits.
(set_slicer_attributes_dvrp)

Ocean atmosphere coupling allows to use independent precursor runs in order
to account for different spin-up times. The time when coupling has to be
started is given by new inipar parameter coupling_start_time. The precursor
ocean run has to be started using new mrun option "-y" in order to add
appendix "_O" to all output files.
(check_for_restart, check_parameters, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, header, init_coupling, modules, mrun,
parin, read_var_list, surface_coupler, time_integration, write_var_list)

Polygon reduction for topography and ground plate isosurface. Reduction level
for buildings can be chosen with parameter cluster_size. (init_dvrp)

External pressure gradient (check_parameters, header, init_3d_model, modules,
parin, prognostic_equations, read_var_list, write_var_list)

New topography case 'single_street_canyon' (header, init_grid, modules, parin,
read_var_list, user_check_parameters, user_header, user_init_grid, write_var_list)

Option to predefine a target bulk velocity for conserve_volume_flow
(check_parameters, header, init_3d_model, modules, parin, read_var_list,
write_var_list)

Option for user defined 2D data output in xy cross sections at z=nzb+1
(data_output_2d, user_data_output_2d)

xy cross section output of surface heatfluxes (latent, sensible)
(average_3d_data, check_parameters, data_output_2d, modules, read_3d_binary,
sum_up_3d_data, write_3d_binary)

average_3d_data, check_for_restart, check_parameters, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, init_coupling, init_dvrp, init_grid, init_3d_model, header, mbuild, modules, mrun, package_parin, parin, prognostic_equations, read_3d_binary, read_var_list, subjob, surface_coupler, timestep, time_integration, user_check_parameters, user_data_output_2d, user_data_output_dvrp, user_header, user_init_grid, write_3d_binary, write_var_list

New: set_particle_attributes, set_slicer_attributes_dvrp

Changed:

lcmuk changed to lc to avoid problems with Intel compiler on sgi-ice
(poisfft)

For extended NetCDF files, the updated title attribute includes an update of
time_average_text where appropriate. (netcdf)

In case of restart runs without extension, initial profiles are not written
to NetCDF-file anymore. (data_output_profiles, modules, read_var_list, write_var_list)

Small change in formatting of the message handling routine concering the output in the
job protocoll. (message)

initializing_actions='read_data_for_recycling' renamed to 'cyclic_fill', now
independent of turbulent_inflow (check_parameters, header, init_3d_model)

2 NetCDF error numbers changed. (data_output_3d)

A Link to the website appendix_a.html is printed for further information
about the possible errors. (message)

Temperature gradient criterion for estimating the boundary layer height
replaced by the gradient criterion of Sullivan et al. (1998). (flow_statistics)

NetCDF unit attribute in timeseries output in case of statistic regions added
(netcdf)

Output of NetCDF messages with aid of message handling routine.
(check_open, close_file, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, netcdf, output_particles_netcdf)

Output of messages replaced by message handling routine.
(advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart,
check_open, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp,
data_output_profiles, data_output_spectra, fft_xy, flow_statistics, header,
init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid,
netcdf, parin, plant_canopy_model, poisfft_hybrid, poismg, read_3d_binary,
read_var_list, surface_coupler, temperton_fft, timestep, user_actions,
user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy,
user_parin, user_read_restart_data, user_spectra )

Maximum number of tails is calculated from maximum number of particles and
skip_particles_for_tail (init_particles)

Value of vertical_particle_advection may differ for each particle group
(advec_particles, header, modules)

First constant in array den also defined as type double. (eqn_state_seawater)

Parameter dvrp_psize moved from particles_par to dvrp_graphics_par. (package_parin)

topography_grid_convention moved from userpar to inipar (check_parameters,
header, parin, read_var_list, user_check_parameters, user_header,
user_init_grid, user_parin, write_var_list)

Default value of grid_matching changed to strict.

Adjustments for runs on lcxt4 (necessary due to an software update on CRAY) and
for coupled runs on ibmy (mrun, subjob)

advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart, check_open, check_parameters, close_file, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, eqn_state_seawater, fft_xy, flow_statistics, header, init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid, message, mrun, netcdf, output_particles_netcdf, package_parin, parin, plant_canopy_model, poisfft, poisfft_hybrid, poismg, read_3d_binary, read_var_list, sort_particles, subjob, user_check_parameters, user_header, user_init_grid, user_parin, surface_coupler, temperton_fft, timestep, user_actions, user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy, user_parin, user_read_restart_data, user_spectra, write_var_list

Errors:

Bugfix: Initial hydrostatic pressure profile in case of ocean runs is now
calculated in 5 iteration steps. (init_ocean)

Bugfix: wrong sign in buoyancy production of ocean part in case of not using
the reference density (only in 3D routine production_e) (production_e)

Bugfix: output of averaged 2d/3d quantities requires that an avaraging
interval has been set, respective error message is included (check_parameters)

Bugfix: Output on unit 14 only if requested by write_binary.
(user_last_actions)

Bugfix to avoid zero division by km_neutral (production_e)

Bugfix for extended NetCDF files: In order to avoid 'data mode' errors if
updated attributes are larger than their original size, NF90_PUT_ATT is called
in 'define mode' enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a
possible performance loss; an alternative strategy would be to ensure equal
attribute size in a job chain. (netcdf)

Bugfix: correction of initial volume flow for non-flat topography (init_3d_model)
Bugfix: zero initialization of arrays within buildings for 'cyclic_fill' (init_3d_model)

Bugfix: to_be_resorted => s_av for time-averaged scalars (data_output_2d, data_output_3d)

Bugfix: error in formatting the output (message)

Bugfix: avoid that ngp_2dh_s_inner becomes zero (init_3_model)

Typographical error: unit of wpt in dots_unit (modules)

Bugfix: error in check, if particles moved further than one subdomain length.
This check must not be applied for newly released particles. (advec_particles)

Bugfix: several tail counters are initialized, particle_tail_coordinates is
only written to file if its third index is > 0, arrays for tails are allocated
with a minimum size of 10 tails if there is no tail initially (init_particles,
advec_particles)

Bugfix: pressure included for profile output (check_parameters)

Bugfix: Type of count and count_rate changed to default INTEGER on NEC machines
(cpu_log)

Bugfix: output if particle time series only if particle advection is switched
on. (time_integration)

Bugfix: qsws was calculated in case of constant heatflux = .FALSE. (prandtl_fluxes)

Bugfix: averaging along z is not allowed for 2d quantities (e.g. u* and z0) (data_output_2d)

Typographical errors (netcdf)

If the inversion height calculated by the prerun is zero, inflow_damping_height
must be explicitly specified (init_3d_model)

Small bugfix concerning 3d 64bit netcdf output format (header)

Bugfix: dt_fixed removed from the restart file, because otherwise, no change
from a fixed to a variable timestep would be possible in restart runs.
(read_var_list, write_var_list)

Bugfix: initial setting of time_coupling in coupled restart runs (time_integration)

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