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

Last change on this file since 132 was 132, checked in by letzel, 14 years ago

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

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

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