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

Last change on this file since 1310 was 1310, checked in by raasch, 7 years ago

update of GPL copyright

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