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

Last change on this file since 3 was 3, checked in by raasch, 15 years ago

RCS Log replace by Id keyword, revision history cleaned up

File size: 4.9 KB
Line 
1 SUBROUTINE diffusivities( theta )
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id$
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                l_stable = 0.76 * sqrt_e(k) / &
93                                  SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
94             ELSE
95                l_stable = l_grid(k)
96             ENDIF
97!
98!--          Adjustment of the mixing length
99             IF ( wall_adjustment )  THEN
100                l(k)  = MIN( l_wall(k,j,i), l_grid(k), l_stable )
101                ll(k) = MIN( l_wall(k,j,i), l_grid(k) )
102             ELSE
103                l(k)  = MIN( l_grid(k), l_stable )
104                ll(k) = l_grid(k)
105             ENDIF
106             IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
107                l(k)  = MIN( l(k),  kappa * zu(k) / phi_m )
108                ll(k) = MIN( ll(k), kappa * zu(k) / phi_m )
109             ENDIF
110
111!
112!--          Compute diffusion coefficients for momentum and heat
113             km(k,j,i) = 0.1 * l(k) * sqrt_e(k)
114             kh(k,j,i) = ( 1.0 + 2.0 * l(k) / ll(k) ) * km(k,j,i)
115
116          ENDDO
117
118!
119!--       Summation for averaged profile (cf. flow_statistics)
120          DO  sr = 0, statistic_regions
121             IF ( rmask(j,i,sr) /= 0.0 )  THEN
122                DO  k = nzb_s_outer(j,i)+1, nzt
123                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l(k)
124                ENDDO
125             ENDIF
126          ENDDO
127
128       ENDDO
129    ENDDO
130
131    sums_l_l(nzt+1,:,tn) = sums_l_l(nzt,:,tn)   ! quasi boundary-condition for
132                                                ! data output
133
134    !$OMP END PARALLEL
135
136!
137!-- Set vertical boundary values (Neumann conditions both at bottom and top).
138!-- Horizontal boundary conditions at vertical walls are not set because
139!-- so far vertical walls require usage of a Prandtl-layer where the boundary
140!-- values of the diffusivities are not needed
141    !$OMP PARALLEL DO
142    DO  i = nxl-1, nxr+1
143       DO  j = nys-1, nyn+1
144          km(nzb_s_inner(j,i),j,i) = km(nzb_s_inner(j,i)+1,j,i)
145          km(nzt+1,j,i)            = km(nzt,j,i)
146          kh(nzb_s_inner(j,i),j,i) = kh(nzb_s_inner(j,i)+1,j,i)
147          kh(nzt+1,j,i)            = kh(nzt,j,i)
148       ENDDO
149    ENDDO
150
151!
152!-- Set Neumann boundary conditions at the outflow boundaries in case of
153!-- non-cyclic lateral boundaries
154    IF ( outflow_l )  THEN
155       km(:,:,nxl-1) = km(:,:,nxl)
156       kh(:,:,nxl-1) = kh(:,:,nxl)
157    ENDIF
158    IF ( outflow_r )  THEN
159       km(:,:,nxr+1) = km(:,:,nxr)
160       kh(:,:,nxr+1) = kh(:,:,nxr)
161    ENDIF
162    IF ( outflow_s )  THEN
163       km(:,nys-1,:) = km(:,nys,:)
164       kh(:,nys-1,:) = kh(:,nys,:)
165    ENDIF
166    IF ( outflow_n )  THEN
167       km(:,nyn+1,:) = km(:,nyn,:)
168       kh(:,nyn+1,:) = kh(:,nyn,:)
169    ENDIF
170
171
172 END SUBROUTINE diffusivities
Note: See TracBrowser for help on using the repository browser.