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

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

Bugfix: summation of sums_l_l in diffusivities.

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