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

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

preliminary uncomplete changes for ocean version

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