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

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

test version of diffusivities.f90 from last revision removed

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