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

Last change on this file since 667 was 667, checked in by suehring, 11 years ago

summary:


Gryschka:

  • Coupling with different resolution and different numbers of PEs in ocean and atmosphere is available
  • Exchange of u and v from ocean surface to atmosphere surface
  • Mirror boundary condition for u and v at the bottom are replaced by dirichlet boundary conditions
  • Inflow turbulence is now defined by flucuations around spanwise mean
  • Bugfixes for cyclic_fill and constant_volume_flow

Suehring:

  • New advection added ( Wicker and Skamarock 5th order ), therefore:
    • New module advec_ws.f90
    • Modified exchange of ghost boundaries.
    • Modified evaluation of turbulent fluxes
    • New index bounds nxlg, nxrg, nysg, nyng

advec_ws.f90


Advection scheme for scalars and momentum using the flux formulation of
Wicker and Skamarock 5th order.
Additionally the module contains of a routine using for initialisation and
steering of the statical evaluation. The computation of turbulent fluxes takes
place inside the advection routines.
In case of vector architectures Dirichlet and Radiation boundary conditions are
outstanding and not available. Furthermore simulations within topography are
not possible so far. A further routine local_diss_ij is available and is used
if a control of dissipative fluxes is desired.

check_parameters.f90


Exchange of parameters between ocean and atmosphere via PE0
Check for illegal combination of ws-scheme and timestep scheme.
Check for topography and ws-scheme.
Check for not cyclic boundary conditions in combination with ws-scheme and
loop_optimization = 'vector'.
Check for call_psolver_at_all_substeps and ws-scheme for momentum_advec.

Different processor/grid topology in atmosphere and ocean is now allowed!
Bugfixes in checking for conserve_volume_flow_mode.

exchange_horiz.f90


Dynamic exchange of ghost points with nbgp_local to ensure that no useless
ghost points exchanged in case of multigrid. type_yz(0) and type_xz(0) used for
normal grid, the remaining types used for the several grid levels.
Exchange is done via MPI-Vectors with a dynamic value of ghost points which
depend on the advection scheme. Exchange of left and right PEs is 10% faster
with MPI-Vectors than without.

flow_statistics.f90


When advection is computed with ws-scheme, turbulent fluxes are already
computed in the respective advection routines and buffered in arrays
sums_xxxx_ws_l(). This is due to a consistent treatment of statistics
with the numerics and to avoid unphysical kinks near the surface. So some if-
requests has to be done to dicern between fluxes from ws-scheme other advection
schemes. Furthermore the computation of z_i is only done if the heat flux
exceeds a minimum value. This affects only simulations of a neutral boundary
layer and is due to reasons of computations in the advection scheme.

inflow_turbulence.f90


Using nbgp recycling planes for a better resolution of the turbulent flow near
the inflow.

init_grid.f90


Definition of new array bounds nxlg, nxrg, nysg, nyng on each PE.
Furthermore the allocation of arrays and steering of loops is done with these
parameters. Call of exchange_horiz are modified.
In case of dirichlet bounday condition at the bottom zu(0)=0.0
dzu_mg has to be set explicitly for a equally spaced grid near bottom.
ddzu_pres added to use a equally spaced grid near bottom.

init_pegrid.f90


Moved determination of target_id's from init_coupling
Determination of parameters needed for coupling (coupling_topology, ngp_a, ngp_o)
with different grid/processor-topology in ocean and atmosphere

Adaption of ngp_xy, ngp_y to a dynamic number of ghost points.
The maximum_grid_level changed from 1 to 0. 0 is the normal grid, 1 to
maximum_grid_level the grids for multigrid, in which 0 and 1 are normal grids.
This distinction is due to reasons of data exchange and performance for the
normal grid and grids in poismg.
The definition of MPI-Vectors adapted to a dynamic numer of ghost points.
New MPI-Vectors for data exchange between left and right boundaries added.
This is due to reasons of performance (10% faster).

ATTENTION: nnz_x undefined problem still has to be solved!!!!!!!!
TEST OUTPUT (TO BE REMOVED) logging mpi2 ierr values

parin.f90


Steering parameter dissipation_control added in inipar.

Makefile


Module advec_ws added.

Modules


Removed u_nzb_p1_for_vfc and v_nzb_p1_for_vfc

For coupling with different resolution in ocean and atmophere:
+nx_a, +nx_o, ny_a, +ny_o, ngp_a, ngp_o, +total_2d_o, +total_2d_a,
+coupling_topology

Buffer arrays for the left sided advective fluxes added in arrays_3d.
+flux_s_u, +flux_s_v, +flux_s_w, +diss_s_u, +diss_s_v, +diss_s_w,
+flux_s_pt, +diss_s_pt, +flux_s_e, +diss_s_e, +flux_s_q, +diss_s_q,
+flux_s_sa, +diss_s_sa
3d arrays for dissipation control added. (only necessary for vector arch.)
+var_x, +var_y, +var_z, +gamma_x, +gamma_y, +gamma_z
Default of momentum_advec and scalar_advec changed to 'ws-scheme' .
+exchange_mg added in control_parameters to steer the data exchange.
Parameters +nbgp, +nxlg, +nxrg, +nysg, +nyng added in indices.
flag array +boundary_flags added in indices to steer the degradation of order
of the advective fluxes when non-cyclic boundaries are used.
MPI-datatypes +type_y, +type_y_int and +type_yz for data_exchange added in
pegrid.
+sums_wsus_ws_l, +sums_wsvs_ws_l, +sums_us2_ws_l, +sums_vs2_ws_l,
+sums_ws2_ws_l, +sums_wspts_ws_l, +sums_wssas_ws_l, +sums_wsqs_ws_l
and +weight_substep added in statistics to steer the statistical evaluation
of turbulent fluxes in the advection routines.
LOGICALS +ws_scheme_sca and +ws_scheme_mom added to get a better performance
in prognostic_equations.
LOGICAL +dissipation_control control added to steer numerical dissipation
in ws-scheme.

Changed length of string run_description_header

pres.f90


New allocation of tend when ws-scheme and multigrid is used. This is due to
reasons of perforance of the data_exchange. The same is done with p after
poismg is called.
nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng when no
multigrid is used. Calls of exchange_horiz are modified.

bugfix: After pressure correction no volume flow correction in case of
non-cyclic boundary conditions
(has to be done only before pressure correction)

Call of SOR routine is referenced with ddzu_pres.

prognostic_equations.f90


Calls of the advection routines with WS5 added.
Calls of ws_statistics added to set the statistical arrays to zero after each
time step.

advec_particles.f90


Declaration of de_dx, de_dy, de_dz adapted to additional ghost points.
Furthermore the calls of exchange_horiz were modified.

asselin_filter.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

average_3d_data.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

boundary_conds.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
Removed mirror boundary conditions for u and v at the bottom in case of
ibc_uv_b == 0. Instead, dirichelt boundary conditions (u=v=0) are set
in init_3d_model

calc_liquid_water_content.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

calc_spectra.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng for
allocation of tend.

check_open.f90


Output of total array size was adapted to nbgp.

data_output_2d.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays local_2d and total_2d.
Calls of exchange_horiz are modified.

data_output_2d.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays. Calls of exchange_horiz are modified.
Skip-value skip_do_avs changed to a dynamic adaption of ghost points.

data_output_mask.f90


Calls of exchange_horiz are modified.

diffusion_e.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_s.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_u.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_v.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_w.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusivities.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusivities.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Calls of exchange_horiz are modified.

exchange_horiz_2d.f90


Dynamic exchange of ghost points with nbgp, which depends on the advection
scheme. Exchange between left and right PEs is now done with MPI-vectors.

global_min_max.f90


Adapting of the index arrays, because MINLOC assumes lowerbound
at 1 and not at nbgp.

init_3d_model.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays. Calls of exchange_horiz are modified.
Call ws_init to initialize arrays needed for statistical evaluation and
optimization when ws-scheme is used.
Initial volume flow is now calculated by using the variable hom_sum.
Therefore the correction of initial volume flow for non-flat topography
removed (removed u_nzb_p1_for_vfc and v_nzb_p1_for_vfc)
Changed surface boundary conditions for u and v in case of ibc_uv_b == 0 from
mirror bc to dirichlet boundary conditions (u=v=0), so that k=nzb is
representative for the height z0

Bugfix: type conversion of '1' to 64bit for the MAX function (ngp_3d_inner)

init_coupling.f90


determination of target_id's moved to init_pegrid

init_pt_anomaly.f90


Call of exchange_horiz are modified.

init_rankine.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Calls of exchange_horiz are modified.

init_slope.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

header.f90


Output of advection scheme.

poismg.f90


Calls of exchange_horiz are modified.

prandtl_fluxes.f90


Changed surface boundary conditions for u and v from mirror bc to dirichelt bc,
therefore u(uzb,:,:) and v(nzb,:,:) is now representative for the height z0
nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

production_e.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

read_3d_binary.f90


+/- 1 replaced with +/- nbgp when swapping and allocating variables.

sor.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Call of exchange_horiz are modified.
bug removed in declaration of ddzw(), nz replaced by nzt+1

subsidence.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

sum_up_3d_data.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

surface_coupler.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in
MPI_SEND() and MPI_RECV.
additional case for nonequivalent processor and grid topopolgy in ocean and
atmosphere added (coupling_topology = 1)

Added exchange of u and v from Ocean to Atmosphere

time_integration.f90


Calls of exchange_horiz are modified.
Adaption to slooping surface.

timestep.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

user_3d_data_averaging.f90, user_data_output_2d.f90, user_data_output_3d.f90,
user_actions.f90, user_init.f90, user_init_plant_canopy.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

user_read_restart_data.f90


Allocation with nbgp.

wall_fluxes.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

write_compressed.f90


Array bounds and nx, ny adapted with nbgp.

sor.f90


bug removed in declaration of ddzw(), nz replaced by nzt+1

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