source: palm/trunk/SOURCE/diffusivities.f90 @ 392

Last change on this file since 392 was 139, checked in by raasch, 14 years ago

New:
---

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

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

For unknown variables (CASE DEFAULT) call new subroutine user_data_output_dvrp

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

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

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

Changed:


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

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

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

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

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

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

non-zero wall_heatflux,

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

humidity/scalar/salinity in a future release.

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

Errors:


Bugfix: summation of sums_l_l in diffusivities.

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

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

Bugfix: model_string needed a default value (combine_plot_fields)

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

Bugfix: assignment of fluxes at walls

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

advec_particles, check_parameters, combine_plot_fields, diffusion_s, diffusivities, init_ocean, init_1d_model, poisfft, production_e

  • Property svn:keywords set to Id
File size: 6.2 KB
Line 
1 SUBROUTINE diffusivities( var, var_reference )
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: diffusivities.f90 139 2007-11-29 09:37:41Z raasch $
11!
12! 137 2007-11-28 08:50:10Z letzel
13! Bugfix for summation of sums_l_l for flow_statistics
14! Vertical scalar profiles now based on nzb_s_inner and ngp_2dh_s_inner.
15!
16! 97 2007-06-21 08:23:15Z raasch
17! Adjustment of mixing length calculation for the ocean version.
18! This is also a bugfix, because the height above the topography is now
19! used instead of the height above level k=0.
20! theta renamed var, dpt_dz renamed dvar_dz, +new argument var_reference
21! use_pt_reference renamed use_reference
22!
23! 57 2007-03-09 12:05:41Z raasch
24! Reference temperature pt_reference can be used in buoyancy term
25!
26! RCS Log replace by Id keyword, revision history cleaned up
27!
28! Revision 1.24  2006/04/26 12:16:26  raasch
29! OpenMP optimization (+sums_l_l_t), sqrt_e must be private
30!
31! Revision 1.1  1997/09/19 07:41:10  raasch
32! Initial revision
33!
34!
35! Description:
36! ------------
37! Computation of the turbulent diffusion coefficients for momentum and heat
38! according to Prandtl-Kolmogorov
39!------------------------------------------------------------------------------!
40
41    USE arrays_3d
42    USE control_parameters
43    USE grid_variables
44    USE indices
45    USE pegrid
46    USE statistics
47
48    IMPLICIT NONE
49
50    INTEGER ::  i, j, k, omp_get_thread_num, sr, tn
51
52    REAL    ::  dvar_dz, l_stable, var_reference
53
54    REAL, SAVE ::  phi_m = 1.0
55
56    REAL    ::  var(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
57
58    REAL, DIMENSION(1:nzt) ::  l, ll, sqrt_e
59
60
61!
62!-- Default thread number in case of one thread
63    tn = 0
64
65!
66!-- Initialization for calculation of the mixing length profile
67    sums_l_l = 0.0
68
69!
70!-- Compute the turbulent diffusion coefficient for momentum
71    !$OMP PARALLEL PRIVATE (dvar_dz,i,j,k,l,ll,l_stable,phi_m,sqrt_e,sr,tn)
72!$  tn = omp_get_thread_num()
73
74    !$OMP DO
75    DO  i = nxl-1, nxr+1
76       DO  j = nys-1, nyn+1
77
78!
79!--       Compute the Phi-function for a possible adaption of the mixing length
80!--       to the Prandtl mixing length
81          IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
82             IF ( rif(j,i) >= 0.0 )  THEN
83                phi_m = 1.0 + 5.0 * rif(j,i)
84             ELSE
85                phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
86             ENDIF
87          ENDIF
88         
89!
90!--       Introduce an optional minimum tke
91          IF ( e_min > 0.0 )  THEN
92             DO  k = nzb_s_inner(j,i)+1, nzt
93                e(k,j,i) = MAX( e(k,j,i), e_min )
94             ENDDO
95          ENDIF
96
97!
98!--       Calculate square root of e in a seperate loop, because it is used
99!--       twice in the next loop (better vectorization)
100          DO  k = nzb_s_inner(j,i)+1, nzt
101             sqrt_e(k) = SQRT( e(k,j,i) )
102          ENDDO
103
104!
105!--       Determine the mixing length
106          DO  k = nzb_s_inner(j,i)+1, nzt
107             dvar_dz = atmos_ocean_sign * &  ! inverse effect of pt/rho gradient
108                       ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
109             IF ( dvar_dz > 0.0 ) THEN
110                IF ( use_reference )  THEN
111                   l_stable = 0.76 * sqrt_e(k) / &
112                                     SQRT( g / var_reference * dvar_dz ) + 1E-5
113                ELSE
114                   l_stable = 0.76 * sqrt_e(k) / &
115                                     SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
116                ENDIF
117             ELSE
118                l_stable = l_grid(k)
119             ENDIF
120!
121!--          Adjustment of the mixing length
122             IF ( wall_adjustment )  THEN
123                l(k)  = MIN( l_wall(k,j,i), l_grid(k), l_stable )
124                ll(k) = MIN( l_wall(k,j,i), l_grid(k) )
125             ELSE
126                l(k)  = MIN( l_grid(k), l_stable )
127                ll(k) = l_grid(k)
128             ENDIF
129             IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
130                l(k)  = MIN( l(k),  kappa * &
131                                    ( zu(k) - zw(nzb_s_inner(j,i)) ) / phi_m )
132                ll(k) = MIN( ll(k), kappa * &
133                                    ( zu(k) - zw(nzb_s_inner(j,i)) ) / phi_m )
134             ENDIF
135
136!
137!--          Compute diffusion coefficients for momentum and heat
138             km(k,j,i) = 0.1 * l(k) * sqrt_e(k)
139             kh(k,j,i) = ( 1.0 + 2.0 * l(k) / ll(k) ) * km(k,j,i)
140
141          ENDDO
142
143!
144!--       Summation for averaged profile (cf. flow_statistics)
145!--       (the IF statement still requires a performance check on NEC machines)
146          DO  sr = 0, statistic_regions
147             IF ( rmask(j,i,sr) /= 0.0 .AND.  &
148                  i >= nxl .AND. i <= nxr .AND. j >= nys .AND. j <= nyn )  THEN
149                DO  k = nzb_s_inner(j,i)+1, nzt
150                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l(k)
151                ENDDO
152             ENDIF
153          ENDDO
154
155       ENDDO
156    ENDDO
157
158    sums_l_l(nzt+1,:,tn) = sums_l_l(nzt,:,tn)   ! quasi boundary-condition for
159                                                ! data output
160
161    !$OMP END PARALLEL
162
163!
164!-- Set vertical boundary values (Neumann conditions both at bottom and top).
165!-- Horizontal boundary conditions at vertical walls are not set because
166!-- so far vertical walls require usage of a Prandtl-layer where the boundary
167!-- values of the diffusivities are not needed
168    !$OMP PARALLEL DO
169    DO  i = nxl-1, nxr+1
170       DO  j = nys-1, nyn+1
171          km(nzb_s_inner(j,i),j,i) = km(nzb_s_inner(j,i)+1,j,i)
172          km(nzt+1,j,i)            = km(nzt,j,i)
173          kh(nzb_s_inner(j,i),j,i) = kh(nzb_s_inner(j,i)+1,j,i)
174          kh(nzt+1,j,i)            = kh(nzt,j,i)
175       ENDDO
176    ENDDO
177
178!
179!-- Set Neumann boundary conditions at the outflow boundaries in case of
180!-- non-cyclic lateral boundaries
181    IF ( outflow_l )  THEN
182       km(:,:,nxl-1) = km(:,:,nxl)
183       kh(:,:,nxl-1) = kh(:,:,nxl)
184    ENDIF
185    IF ( outflow_r )  THEN
186       km(:,:,nxr+1) = km(:,:,nxr)
187       kh(:,:,nxr+1) = kh(:,:,nxr)
188    ENDIF
189    IF ( outflow_s )  THEN
190       km(:,nys-1,:) = km(:,nys,:)
191       kh(:,nys-1,:) = kh(:,nys,:)
192    ENDIF
193    IF ( outflow_n )  THEN
194       km(:,nyn+1,:) = km(:,nyn,:)
195       kh(:,nyn+1,:) = kh(:,nyn,:)
196    ENDIF
197
198
199 END SUBROUTINE diffusivities
Note: See TracBrowser for help on using the repository browser.