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

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

preliminary update of further changes, advec_particles is not running!

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