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

Last change on this file since 509 was 509, checked in by raasch, 12 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 6.4 KB
Line 
1 SUBROUTINE diffusivities( var, var_reference )
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: diffusivities.f90 509 2010-03-12 00:31:47Z raasch $
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                   IF ( var(k,j,i) <= 0.0 )  THEN
115                      WRITE (9,*) 'i=', i, ' j=', j, ' k=', k, ' var=', &
116                                  var(k,j,i)
117                      CALL local_flush( 9 )
118                   ENDIF
119                   l_stable = 0.76 * sqrt_e(k) / &
120                                     SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
121                ENDIF
122             ELSE
123                l_stable = l_grid(k)
124             ENDIF
125!
126!--          Adjustment of the mixing length
127             IF ( wall_adjustment )  THEN
128                l(k)  = MIN( l_wall(k,j,i), l_grid(k), l_stable )
129                ll(k) = MIN( l_wall(k,j,i), l_grid(k) )
130             ELSE
131                l(k)  = MIN( l_grid(k), l_stable )
132                ll(k) = l_grid(k)
133             ENDIF
134             IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
135                l(k)  = MIN( l(k),  kappa * &
136                                    ( zu(k) - zw(nzb_s_inner(j,i)) ) / phi_m )
137                ll(k) = MIN( ll(k), kappa * &
138                                    ( zu(k) - zw(nzb_s_inner(j,i)) ) / phi_m )
139             ENDIF
140
141!
142!--          Compute diffusion coefficients for momentum and heat
143             km(k,j,i) = 0.1 * l(k) * sqrt_e(k)
144             kh(k,j,i) = ( 1.0 + 2.0 * l(k) / ll(k) ) * km(k,j,i)
145
146          ENDDO
147
148!
149!--       Summation for averaged profile (cf. flow_statistics)
150!--       (the IF statement still requires a performance check on NEC machines)
151          DO  sr = 0, statistic_regions
152             IF ( rmask(j,i,sr) /= 0.0 .AND.  &
153                  i >= nxl .AND. i <= nxr .AND. j >= nys .AND. j <= nyn )  THEN
154                DO  k = nzb_s_inner(j,i)+1, nzt
155                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l(k)
156                ENDDO
157             ENDIF
158          ENDDO
159
160       ENDDO
161    ENDDO
162
163    sums_l_l(nzt+1,:,tn) = sums_l_l(nzt,:,tn)   ! quasi boundary-condition for
164                                                ! data output
165
166    !$OMP END PARALLEL
167
168!
169!-- Set vertical boundary values (Neumann conditions both at bottom and top).
170!-- Horizontal boundary conditions at vertical walls are not set because
171!-- so far vertical walls require usage of a Prandtl-layer where the boundary
172!-- values of the diffusivities are not needed
173    !$OMP PARALLEL DO
174    DO  i = nxl-1, nxr+1
175       DO  j = nys-1, nyn+1
176          km(nzb_s_inner(j,i),j,i) = km(nzb_s_inner(j,i)+1,j,i)
177          km(nzt+1,j,i)            = km(nzt,j,i)
178          kh(nzb_s_inner(j,i),j,i) = kh(nzb_s_inner(j,i)+1,j,i)
179          kh(nzt+1,j,i)            = kh(nzt,j,i)
180       ENDDO
181    ENDDO
182
183!
184!-- Set Neumann boundary conditions at the outflow boundaries in case of
185!-- non-cyclic lateral boundaries
186    IF ( outflow_l )  THEN
187       km(:,:,nxl-1) = km(:,:,nxl)
188       kh(:,:,nxl-1) = kh(:,:,nxl)
189    ENDIF
190    IF ( outflow_r )  THEN
191       km(:,:,nxr+1) = km(:,:,nxr)
192       kh(:,:,nxr+1) = kh(:,:,nxr)
193    ENDIF
194    IF ( outflow_s )  THEN
195       km(:,nys-1,:) = km(:,nys,:)
196       kh(:,nys-1,:) = kh(:,nys,:)
197    ENDIF
198    IF ( outflow_n )  THEN
199       km(:,nyn+1,:) = km(:,nyn,:)
200       kh(:,nyn+1,:) = kh(:,nyn,:)
201    ENDIF
202
203
204 END SUBROUTINE diffusivities
Note: See TracBrowser for help on using the repository browser.