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

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

ONLY-attribute added to USE-statements,
kind-parameters added to all INTEGER and REAL declaration statements,
kinds are defined in new module kinds,
old module precision_kind is removed,
revision history before 2012 removed,
comment fields (!:) to be used for variable explanations added to all variable declaration statements

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