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

Last change on this file since 128 was 98, checked in by raasch, 17 years ago

updating comments and rc-file

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