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

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

Initial repository layout and content

File size: 7.9 KB
Line 
1 SUBROUTINE diffusivities( theta )
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Log: diffusivities.f90,v $
11! Revision 1.24  2006/04/26 12:16:26  raasch
12! OpenMP optimization (+sums_l_l_t), sqrt_e must be private
13!
14! Revision 1.23  2006/02/23 12:13:29  raasch
15! nzb_2d replaced by nzb_s_inner/outer, option of a minimum tke
16! quasi boundary conditions is now set for sums_l_l instead of sums_l
17!
18! Revision 1.22  2005/03/26 20:16:26  raasch
19! Eddy diffusivities are also calculated for the ghost points in order to spare
20! the communication otherwise needed for the exchange of these points.
21! Neumann boundary conditions set at the outflow boundaries in case of
22! non-cyclic lateral boundaries
23!
24! Revision 1.21  2004/01/30 10:23:12  raasch
25! Scalar lower k index nzb replaced by 2d-array nzb_2d
26!
27! Revision 1.20  2003/03/14 13:40:09  raasch
28! Loop optimization (SQRT(e) solved in a seperate loop, joining of loops)
29!
30! Revision 1.19  2003/03/12 16:28:56  raasch
31! phi_m is initialized in declaration statement in order to avoid run time
32! errors due to compiler problems on ibm and nec
33!
34! Revision 1.18  2002/12/19 14:32:23  raasch
35! Correction of mixing length term (l(k)/ll(k)). The condition kh=3*km in
36! the unstable case is now also exactly met in the wall adjustment region.
37! Factor 0.7 in wall adjustment part replaced by variable
38! wall_adjustment_factor, which is set to 1.8 in modules.f90.
39! Missing OMP DO statement added.
40!
41! Revision 1.17  2002/06/11 12:56:49  raasch
42! OpenMP directives added,
43! array l changed from allocatable to automatic, because it must be on the stack
44! in case of using OpenMP
45!
46! Revision 1.16  2001/08/21  08:30:34  08:30:34  raasch (Siegfried Raasch)
47! Wall adjustment of mixing length to 0.7 z can be switched off
48! Lower boundary condition for km and kh changed to Neumann
49!
50! Revision 1.15  2001/03/30 07:22:09  raasch
51! Near surface mixing length is limited to 0.7*zu,
52! if adjust_mixing_length=T, l is modified at all gridpoints (previously at
53! k=1 only),
54! Translation of remaining German identifiers (variables, subroutines, etc.)
55!
56! Revision 1.14  2001/01/22 06:31:36  raasch
57! Module test_variables removed
58!
59! Revision 1.13  2001/01/02 17:27:41  raasch
60! -dpt_dz_d, dpt_dz_u
61!
62! Revision 1.12  2000/12/20 10:35:47  letzel
63! All comments translated into English.
64!
65! Revision 1.11  2000/04/18 08:10:42  schroeter
66! Revision 1.4 wieder rueckgaengig gemacht, das Stabilitaets-
67! kriterium basiert nun wieder auf zentralen Differenzen
68!
69! Revision 1.9  2000/04/13 13:42:14  schroeter
70! je nach Initialisierungsmodus (trocken/feucht) fliesst in
71! Mischungswegberechnung pt oder vpt ein, wird durch entsprechende
72! Variablenuebergabe geregelt
73!
74! Revision 1.8  99/02/17  09:17:39  09:17:39  raasch (Siegfried Raasch)
75! Kriterium fuer reduzierten Mischungsweg im stabil geschichteten Fall enger
76! gefasst (Schichtung muss sowohl oberhalb als auch unterhalb des betrachteten
77! Gitterpunkts stabil sein)
78!
79! Revision 1.7  1998/09/22 17:18:41  raasch
80! Testweise cm = 0.4 bei k=1, aber vorerst auskommentiert
81!
82! Revision 1.6  1998/08/05 06:52:44  raasch
83! Mischungsweganpassung an Prandtl nur bei k=1
84!
85! Revision 1.2  1998/03/11 11:49:26  raasch
86! Anpassung des Mischungsweges an den Prandtlschen Mischungsweg moeglich
87!
88! Revision 1.1  1997/09/19 07:41:10  raasch
89! Initial revision
90!
91!
92! Description:
93! ------------
94! Computation of the turbulent diffusion coefficients for momentum and heat
95! according to Prandtl-Kolmogorov
96!------------------------------------------------------------------------------!
97
98    USE arrays_3d
99    USE control_parameters
100    USE grid_variables
101    USE indices
102    USE pegrid
103    USE statistics
104
105    IMPLICIT NONE
106
107    INTEGER ::  i, j, k, omp_get_thread_num, sr, tn
108
109    REAL    ::  dpt_dz, l_stable, phi_m = 1.0
110
111    REAL    ::  theta(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
112
113    REAL, DIMENSION(1:nzt) ::  l, ll, sqrt_e
114
115
116!
117!-- Default thread number in case of one thread
118    tn = 0
119
120!
121!-- Initialization for calculation of the mixing length profile
122    sums_l_l = 0.0
123
124!
125!-- Compute the turbulent diffusion coefficient for momentum
126    !$OMP PARALLEL PRIVATE (dpt_dz,i,j,k,l,ll,l_stable,phi_m,sqrt_e,sr,tn)
127!$  tn = omp_get_thread_num()
128
129    !$OMP DO
130    DO  i = nxl-1, nxr+1
131       DO  j = nys-1, nyn+1
132
133!
134!--       Compute the Phi-function for a possible adaption of the mixing length
135!--       to the Prandtl mixing length
136          IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
137             IF ( rif(j,i) >= 0.0 )  THEN
138                phi_m = 1.0 + 5.0 * rif(j,i)
139             ELSE
140                phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
141             ENDIF
142          ENDIF
143         
144!
145!--       Introduce an optional minimum tke
146          IF ( e_min > 0.0 )  THEN
147             DO  k = nzb_s_inner(j,i)+1, nzt
148                e(k,j,i) = MAX( e(k,j,i), e_min )
149             ENDDO
150          ENDIF
151
152!
153!--       Calculate square root of e in a seperate loop, because it is used
154!--       twice in the next loop (better vectorization)
155          DO  k = nzb_s_inner(j,i)+1, nzt
156             sqrt_e(k) = SQRT( e(k,j,i) )
157          ENDDO
158
159!
160!--       Determine the mixing length
161          DO  k = nzb_s_inner(j,i)+1, nzt
162             dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
163             IF ( dpt_dz > 0.0 ) THEN
164                l_stable = 0.76 * sqrt_e(k) / &
165                                  SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
166             ELSE
167                l_stable = l_grid(k)
168             ENDIF
169!
170!--          Adjustment of the mixing length
171             IF ( wall_adjustment )  THEN
172                l(k)  = MIN( l_wall(k,j,i), l_grid(k), l_stable )
173                ll(k) = MIN( l_wall(k,j,i), l_grid(k) )
174             ELSE
175                l(k)  = MIN( l_grid(k), l_stable )
176                ll(k) = l_grid(k)
177             ENDIF
178             IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
179                l(k)  = MIN( l(k),  kappa * zu(k) / phi_m )
180                ll(k) = MIN( ll(k), kappa * zu(k) / phi_m )
181             ENDIF
182
183!
184!--          Compute diffusion coefficients for momentum and heat
185             km(k,j,i) = 0.1 * l(k) * sqrt_e(k)
186             kh(k,j,i) = ( 1.0 + 2.0 * l(k) / ll(k) ) * km(k,j,i)
187
188          ENDDO
189
190!
191!--       Summation for averaged profile (cf. flow_statistics)
192          DO  sr = 0, statistic_regions
193             IF ( rmask(j,i,sr) /= 0.0 )  THEN
194                DO  k = nzb_s_outer(j,i)+1, nzt
195                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l(k)
196                ENDDO
197             ENDIF
198          ENDDO
199
200       ENDDO
201    ENDDO
202
203    sums_l_l(nzt+1,:,tn) = sums_l_l(nzt,:,tn)   ! quasi boundary-condition for
204                                                ! data output
205
206    !$OMP END PARALLEL
207
208!
209!-- Set vertical boundary values (Neumann conditions both at bottom and top).
210!-- Horizontal boundary conditions at vertical walls are not set because
211!-- so far vertical walls require usage of a Prandtl-layer where the boundary
212!-- values of the diffusivities are not needed
213    !$OMP PARALLEL DO
214    DO  i = nxl-1, nxr+1
215       DO  j = nys-1, nyn+1
216          km(nzb_s_inner(j,i),j,i) = km(nzb_s_inner(j,i)+1,j,i)
217          km(nzt+1,j,i)            = km(nzt,j,i)
218          kh(nzb_s_inner(j,i),j,i) = kh(nzb_s_inner(j,i)+1,j,i)
219          kh(nzt+1,j,i)            = kh(nzt,j,i)
220       ENDDO
221    ENDDO
222
223!
224!-- Set Neumann boundary conditions at the outflow boundaries in case of
225!-- non-cyclic lateral boundaries
226    IF ( outflow_l )  THEN
227       km(:,:,nxl-1) = km(:,:,nxl)
228       kh(:,:,nxl-1) = kh(:,:,nxl)
229    ENDIF
230    IF ( outflow_r )  THEN
231       km(:,:,nxr+1) = km(:,:,nxr)
232       kh(:,:,nxr+1) = kh(:,:,nxr)
233    ENDIF
234    IF ( outflow_s )  THEN
235       km(:,nys-1,:) = km(:,nys,:)
236       kh(:,nys-1,:) = kh(:,nys,:)
237    ENDIF
238    IF ( outflow_n )  THEN
239       km(:,nyn+1,:) = km(:,nyn,:)
240       kh(:,nyn+1,:) = kh(:,nyn,:)
241    ENDIF
242
243
244 END SUBROUTINE diffusivities
Note: See TracBrowser for help on using the repository browser.