Home


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

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

Last change on this file since 91 was 77, checked in by raasch, 18 years ago

New:
---

particle reflection from vertical walls implemented, particle SGS model adjusted to walls

Wall functions for vertical walls now include diabatic conditions. New subroutines wall_fluxes, wall_fluxes_e. New 4D-array rif_wall.

new d3par-parameter netcdf_64bit_3d to switch on 64bit offset only for 3D files

new d3par-parameter dt_max to define the maximum value for the allowed timestep

new inipar-parameter loop_optimization to control the loop optimization method

new inipar-parameter pt_refrence. If given, this value is used as the reference that in buoyancy terms (otherwise, the instantaneous horizontally averaged temperature is used).

new user interface user_advec_particles

new initializing action "by_user" calls user_init_3d_model and allows the initial setting of all 3d arrays

topography height informations are stored on arrays zu_s_inner and zw_w_inner and output to the 2d/3d NetCDF files

samples added to the user interface which show how to add user-define time series quantities.

calculation/output of precipitation amount, precipitation rate and z0 (by setting "pra*", "prr*", "z0*" with data_output). The time interval on which the precipitation amount is defined is set by new d3par-parameter precipitation_amount_interval

unit 9 opened for debug output (file DEBUG_<pe#>)

Makefile, advec_particles, average_3d_data, buoyancy, calc_precipitation, check_open, check_parameters, data_output_2d, diffusion_e, diffusion_u, diffusion_v, diffusion_w, diffusivities, header, impact_of_latent_heat, init_particles, init_3d_model, modules, netcdf, parin, production_e, read_var_list, read_3d_binary, sum_up_3d_data, user_interface, write_var_list, write_3d_binary

New: wall_fluxes

Changed:

General revision of non-cyclic horizontal boundary conditions: radiation boundary conditions are now used instead of Neumann conditions at the outflow (calculation needs velocity values for t-dt, which are stored on new arrays u_m_l, u_m_r, etc.), calculation of mean outflow is not needed any more, volume flow control is added for the outflow boundary (currently only for the north boundary!!), additional gridpoints along x and y (uxrp, vynp) are not needed any more, routine "boundary_conds" now operates on timelevel t+dt and is not split in two parts (main, uvw_outflow) any more, Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary conditions for all 2d-arrays that are handled by exchange_horiz_2d

The FFT-method for solving the Poisson-equation is now working with Neumann boundary conditions both at the bottom and the top. This requires adjustments of the tridiagonal coefficients and subtracting the horizontally averaged mean from the vertical velocity field.

+age_m in particle_type

Particles-package is now part of the default code ("-p particles" is not needed any more).

Move call of user_actions( 'after_integration' ) below increment of times
and counters. user_actions is now called for each statistic region and has as an argument the number of the respective region (sr)

d3par-parameter data_output_ts removed. Timeseries output for "profil" removed. Timeseries are now switched on by dt_dots. Timeseries data is collected in flow_statistics.

Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc.

q is not allowed to become negative (prognostic_equations).

poisfft_init is only called if fft-solver is switched on (init_pegrid).

d3par-parameter moisture renamed to humidity.

Subversion global revision number is read from mrun and added to the run description header and to the run control (_rc) file.

vtk directives removed from main program.

The uitility routine interpret_config reads PALM environment variables from NAMELIST instead using the system call GETENV.

advec_u_pw, advec_u_up, advec_v_pw, advec_v_up, asselin_filter, check_parameters, coriolis, data_output_dvrp, data_output_ptseries, data_output_ts, data_output_2d, data_output_3d, diffusion_u, diffusion_v, exchange_horiz, exchange_horiz_2d, flow_statistics, header, init_grid, init_particles, init_pegrid, init_rankine, init_pt_anomaly, init_1d_model, init_3d_model, modules, palm, package_parin, parin, poisfft, poismg, prandtl_fluxes, pres, production_e, prognostic_equations, read_var_list, read_3d_binary, sor, swap_timelevel, time_integration, write_var_list, write_3d_binary

Errors:

Bugfix: preset of tendencies te_em, te_um, te_vm in init_1d_model

Bugfix in sample for reading user defined data from restart file (user_init)

Bugfix in setting diffusivities for cases with the outflow damping layer extending over more than one subdomain (init_3d_model)

Check for possible negative humidities in the initial humidity profile.

in Makefile, default suffixes removed from the suffix list to avoid calling of m2c in
# case of .mod files

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