source: palm/tags/release-3.2b/SOURCE/diffusivities.f90 @ 4003

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