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

Last change on this file since 1179 was 1179, checked in by raasch, 8 years ago

New:
---
Initial profiles can be used as reference state in the buoyancy term. New parameter
reference_state introduced. Calculation and handling of reference state in buoyancy term revised.
binary version for restart files changed from 3.9 to 3.9a (no downward compatibility!),
initial profile for rho added to hom (id=77)

Errors:


small bugfix for background communication (time_integration)

  • Property svn:keywords set to Id
File size: 6.8 KB
Line 
1 SUBROUTINE diffusivities( var, var_reference )
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2012  Leibniz University Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22! use_reference renamed use_single_reference_value
23!
24! Former revisions:
25! -----------------
26! $Id: diffusivities.f90 1179 2013-06-14 05:57:58Z raasch $
27!
28! 1036 2012-10-22 13:43:42Z raasch
29! code put under GPL (PALM 3.9)
30!
31! 1015 2012-09-27 09:23:24Z raasch
32! OpenACC statements added + code changes required for GPU optimization,
33! adjustment of mixing length to the Prandtl mixing length at first grid point
34! above ground removed
35!
36! 667 2010-12-23 12:06:00Z suehring/gryschka
37! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
38!
39! 137 2007-11-28 08:50:10Z letzel
40! Bugfix for summation of sums_l_l for flow_statistics
41! Vertical scalar profiles now based on nzb_s_inner and ngp_2dh_s_inner.
42!
43! 97 2007-06-21 08:23:15Z raasch
44! Adjustment of mixing length calculation for the ocean version.
45! This is also a bugfix, because the height above the topography is now
46! used instead of the height above level k=0.
47! theta renamed var, dpt_dz renamed dvar_dz, +new argument var_reference
48! use_pt_reference renamed use_reference
49!
50! 57 2007-03-09 12:05:41Z raasch
51! Reference temperature pt_reference can be used in buoyancy term
52!
53! RCS Log replace by Id keyword, revision history cleaned up
54!
55! Revision 1.24  2006/04/26 12:16:26  raasch
56! OpenMP optimization (+sums_l_l_t), sqrt_e must be private
57!
58! Revision 1.1  1997/09/19 07:41:10  raasch
59! Initial revision
60!
61!
62! Description:
63! ------------
64! Computation of the turbulent diffusion coefficients for momentum and heat
65! according to Prandtl-Kolmogorov
66!------------------------------------------------------------------------------!
67
68    USE arrays_3d
69    USE control_parameters
70    USE grid_variables
71    USE indices
72    USE pegrid
73    USE statistics
74
75    IMPLICIT NONE
76
77    INTEGER ::  i, j, k, omp_get_thread_num, sr, tn
78
79    REAL    ::  dvar_dz, l, ll, l_stable, sqrt_e, var_reference
80
81    REAL    ::  var(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
82
83
84!
85!-- Default thread number in case of one thread
86    tn = 0
87
88!
89!-- Initialization for calculation of the mixing length profile
90    sums_l_l = 0.0
91
92!
93!-- Compute the turbulent diffusion coefficient for momentum
94    !$OMP PARALLEL PRIVATE (dvar_dz,i,j,k,l,ll,l_stable,sqrt_e,sr,tn)
95!$  tn = omp_get_thread_num()
96
97!
98!-- Data declerations for accelerators
99    !$acc data present( dd2zu, e, km, kh, l_grid, l_wall, nzb_s_inner, rif, var )
100    !$acc kernels
101
102!
103!-- Introduce an optional minimum tke
104    IF ( e_min > 0.0 )  THEN
105       !$OMP DO
106       !$acc loop
107       DO  i = nxlg, nxrg
108          DO  j = nysg, nyng
109             !$acc loop vector( 32 )
110             DO  k = 1, nzt
111                IF ( k > nzb_s_inner(j,i) )  THEN
112                   e(k,j,i) = MAX( e(k,j,i), e_min )
113                ENDIF
114             ENDDO
115          ENDDO
116       ENDDO
117    ENDIF
118
119    !$OMP DO
120    !$acc loop
121    DO  i = nxlg, nxrg
122       DO  j = nysg, nyng
123          !$acc loop vector( 32 )
124          DO  k = 1, nzt
125
126             IF ( k > nzb_s_inner(j,i) )  THEN
127
128                sqrt_e = SQRT( e(k,j,i) )
129!
130!--             Determine the mixing length
131                dvar_dz = atmos_ocean_sign * &  ! inverse effect of pt/rho gradient
132                          ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
133                IF ( dvar_dz > 0.0 ) THEN
134                   IF ( use_single_reference_value )  THEN
135                      l_stable = 0.76 * sqrt_e / &
136                                 SQRT( g / var_reference * dvar_dz ) + 1E-5
137                   ELSE
138                      l_stable = 0.76 * sqrt_e / &
139                                 SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
140                   ENDIF
141                ELSE
142                   l_stable = l_grid(k)
143                ENDIF
144!
145!--             Adjustment of the mixing length
146                IF ( wall_adjustment )  THEN
147                   = MIN( l_wall(k,j,i), l_grid(k), l_stable )
148                   ll = MIN( l_wall(k,j,i), l_grid(k) )
149                ELSE
150                   = MIN( l_grid(k), l_stable )
151                   ll = l_grid(k)
152                ENDIF
153
154      !
155      !--       Compute diffusion coefficients for momentum and heat
156                km(k,j,i) = 0.1 * l * sqrt_e
157                kh(k,j,i) = ( 1.0 + 2.0 * l / ll ) * km(k,j,i)
158
159             ENDIF
160
161#if ! defined( __openacc )
162!
163!++          Statistics still have to be realized for accelerators
164!--          Summation for averaged profile (cf. flow_statistics)
165             DO  sr = 0, statistic_regions
166                sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l * rmask(j,i,sr)
167             ENDDO
168#endif
169
170          ENDDO
171       ENDDO
172    ENDDO
173
174#if ! defined( __openacc )
175!
176!++ Statistics still have to be realized for accelerators
177    sums_l_l(nzt+1,:,tn) = sums_l_l(nzt,:,tn)   ! quasi boundary-condition for
178                                                  ! data output
179#endif
180    !$OMP END PARALLEL
181
182!
183!-- Set vertical boundary values (Neumann conditions both at bottom and top).
184!-- Horizontal boundary conditions at vertical walls are not set because
185!-- so far vertical walls require usage of a Prandtl-layer where the boundary
186!-- values of the diffusivities are not needed
187    !$OMP PARALLEL DO
188    !$acc loop
189    DO  i = nxlg, nxrg
190       DO  j = nysg, nyng
191          km(nzb_s_inner(j,i),j,i) = km(nzb_s_inner(j,i)+1,j,i)
192          km(nzt+1,j,i)            = km(nzt,j,i)
193          kh(nzb_s_inner(j,i),j,i) = kh(nzb_s_inner(j,i)+1,j,i)
194          kh(nzt+1,j,i)            = kh(nzt,j,i)
195       ENDDO
196    ENDDO
197
198!
199!-- Set Neumann boundary conditions at the outflow boundaries in case of
200!-- non-cyclic lateral boundaries
201    IF ( outflow_l )  THEN
202       km(:,:,nxl-1) = km(:,:,nxl)
203       kh(:,:,nxl-1) = kh(:,:,nxl)
204    ENDIF
205    IF ( outflow_r )  THEN
206       km(:,:,nxr+1) = km(:,:,nxr)
207       kh(:,:,nxr+1) = kh(:,:,nxr)
208    ENDIF
209    IF ( outflow_s )  THEN
210       km(:,nys-1,:) = km(:,nys,:)
211       kh(:,nys-1,:) = kh(:,nys,:)
212    ENDIF
213    IF ( outflow_n )  THEN
214       km(:,nyn+1,:) = km(:,nyn,:)
215       kh(:,nyn+1,:) = kh(:,nyn,:)
216    ENDIF
217
218    !$acc end kernels
219    !$acc end data
220
221 END SUBROUTINE diffusivities
Note: See TracBrowser for help on using the repository browser.